source("Welfare_auxfunctions_simplespouse.R")
source("psoptim_parallel3.R")


transfers<- function(simdata, kinks, hours=2000,  agelength = 1, 
                fsgross = 23920,fsnet = 18400, SSIval = 5100,  FSgenerosity = 1, SSgenerosity = 1) {
  simdata[,SSgenerosity:=SSgenerosity]
  simdata[time >= 40 * agelength,SSgenerosity:=1]   #SSgenerosity only affects DI receipt--not retirement.
  simdata[,sswage:=Wage*hours*1/exp(wagepars$wagereg[[Type]]$health[Health]),by=Type]
  simdata[,sswagesp:=spWage*hours*1/exp(spwagepars$wagereg[[Type]]$health[spHealth]),by=Type]
  simdata[(sswage >= 62700), sswage := 62700]
  simdata[(sswagesp >= 62700), sswagesp := 62700]
  

  simdata[,SSamount:=0]
  simdata[yearDI==1,SSamount:= pmax(0.9 * pmin(kinks[1], sswage) +
                          0.32 * pmin(pmax(0.0, sswage - kinks[1]), kinks[2] - kinks[1]) +
                          0.15 * pmin(pmax(0.0, sswage - kinks[2]), kinks[3] - kinks[2]), SSIval)]
  simdata[,SSamount := SSamount * SSgenerosity]
  simdata[,SStaxable := 0.0]
      simdata[single==0&SSamount>0&(SSamount * 0.5  <= 44000), SStaxable := pmax((SSamount * 0.5 - 32000) * 0.5, 0.0)]
      simdata[single==0&SSamount>0&(SSamount * 0.5 > 44000), SStaxable := (12000) * 0.5]
      simdata[single==0&SSamount>0& (SSamount * 0.5 < SStaxable), SStaxable := SSamount * 0.5]
      simdata[single==0&SSamount>0& (SSamount * 0.5  > 44000), SStaxable := SStaxable + 0.85 * (SSamount * 0.5 - 44000)]
      simdata[single==0&SSamount>0& (SSamount * 0.85 < SStaxable), SStaxable := SSamount * 0.85]

      simdata[single==1&SSamount>0&(SSamount * 0.5  <= 34000), SStaxable := pmax((SSamount * 0.5 - 25000) * 0.5, 0.0)]
      simdata[single==1&SSamount>0&(SSamount * 0.5  > 34000), SStaxable := (9000) * 0.5]
      simdata[single==1&SSamount>0&(SSamount * 0.5 < SStaxable), SStaxable := SSamount * 0.5]
      simdata[single==1&SSamount>0&(SSamount * 0.5  > 34000), SStaxable := SStaxable + 0.85 * (SSamount * 0.5 + pretaxinc - 34000)]
      simdata[single==1&SSamount>0&(SSamount * 0.85 < SStaxable), SStaxable := SSamount * 0.85]

      simdata[,taxableinc := pmax(SStaxable - 6700 * (1 - single) - 4000 * single - 2550 * (1 + (1 - single) + (1-single)*nkids), 0.0)]
  
    simdata[single==0& taxableinc <= 40100, taxsim := pmax(taxableinc*0.15,0)]
    simdata[single==0& taxableinc > 40100 & taxableinc<=96900, taxsim := 40100 * 0.15 + (taxableinc - 40100) * 0.28]
    simdata[single==0&taxableinc>96900 & taxableinc<=147700, taxsim := 40100 * 0.15 + (96900 - 40100) * 0.28 + (taxableinc - 96900) * 0.31]
    simdata[single==0&taxableinc>147700&taxableinc<=263750, taxsim:=40100 * 0.15 + (96900 - 40100) * 0.28 + (147700 - 96900) * 0.31 + (taxableinc - 147700) * 0.36]
    simdata[single==0&taxableinc>263750, taxsim:=40100 * 0.15 + (96900 - 40100) * 0.28 + (147700 - 96900) * 0.31 + (263750 - 147700) * 0.36 + (taxableinc - 263750) * 0.396]
    
    simdata[single==1&(taxableinc <= 24000), taxsim := pmax(taxableinc * 0.15, 0.0)]        
    simdata[single==1&(taxableinc > 24000 & taxableinc <= 58150), taxsim := 24000 * 0.15 + (taxableinc - 24000) * 0.28]
    simdata[single==1&(taxableinc > 58150 & taxableinc <= 121300), taxsim := 24000 * 0.15 + (58150 - 24000) * 0.28 + (taxableinc - 58150) * 0.31]
    simdata[single==1& (taxableinc > 121300 & taxableinc <= 263750), taxsim := 24000 * 0.15 + (58150 - 24000) * 0.28 + (121300 - 58150) * 0.31 + (taxableinc - 121300) * 0.36]
    simdata[single==1& (taxableinc > 263750), taxsim := 24000 * 0.15 + (58150 - 24000) * 0.28 + (121300 - 58150) * 0.31 + (263750 - 121300) * 0.36 + (taxableinc - 263750) * 0.396]
  #adding child tax credit:
    simdata[,taxsim := pmax(taxsim - 400 * nkids * (1-single), 0.0)]
  
  #Including FICA taxes: NONE, because no work

  #INCLUDING EITC: NONE, because no work

  #including food stamps:
    simdata[,potential_transfer:= pmax(0,potential_foodstamps*FSgenerosity - SSamount*0.3)+SSamount]
}


drawselect<-function(simnum){
  draws<-list()
  for(ty in 1:length(simnum)){
    drawnums<-sample(1:nrow(obsdraws[type_fix==ty,]),size=simnum[ty],replace=TRUE)
    draws[[ty]]<-as.vector(t(obsdraws[type_fix==ty,][drawnums,-1]))
  }
  return(unlist(draws))
}

drawselectC<-function(simnum){
  draws<-list()
  for(ty in 1:length(simnum)){
    drawnums<-sample(1:nrow(obsdraws_cons[type_fix==ty,]),size=simnum[ty],replace=TRUE)
    draws[[ty]]<-as.vector(t(obsdraws_cons[type_fix==ty,][drawnums,-1]))
  }
  return(unlist(draws))
}
#Previous values:
#NO RESTRICTIONS on head or spouses in sample, aside from those from LP:
#FEs come specifically from married heads
#Imposing that missing spouses have 10th percentile FE
wagesetup_fit <- function(pars,shift=0,spshift=0,wage_gridsize=10, spwage_gridsize=10,
                      wagereg_agemult=1,
                      spwagereg_agemult=1,covs=TRUE,post=FALSE){
  wagevals<-as.matrix(pars[c("mod","sev","age","age2","Type1","Type2","Type3")])
  spwagevals<-as.matrix(pars[c("spmod","spage","spage2","spType1","spType2","spType3")])
  rownames(spwagevals)<-gsub("sp","",rownames(spwagevals))
  numtypes<-3
  
  wagevars<-as.matrix(c(pars[c("wagevar","spwagevar","wagecorr","noisevar","spnoisevar")]))
  
  if(post==FALSE & covs == TRUE){
    wagepars <- list(wagevar = wagevars[1,1]^0.5/(agelength^0.5),
                     wagecorr = wagevars[3,1],
                     noisevar= wagevars[4,1]^0.5/(agelength^0.5)
    )
    wagepars <-c(wagepars,
                 wagereg = list(lapply(1:numtypes,
                                       function(x) list(cons = 0+shift*wagepars$wagevar+wagevals[paste0("Type",x),1], age = wagevals["age",1]*wagereg_agemult,agesq=wagevals["age2",1]/100, health=c(0,wagevals["mod",1],wagevals["sev",1]))
                 )),
                 agelength = agelength
    )
    
    spwagepars <- list(wagevar = wagevars[2,1]^0.5/(agelength^0.5),
                       wagecorr = wagevars[3,1],
                       noisevar= wagevars[5,1]^0.5/(agelength^0.5)
    )
    spwagepars<-c(spwagepars,
                  wagereg = list(lapply(1:numtypes,
                                        function(x) list(cons = 0+spshift*spwagepars$wagevar+spwagevals[paste0("Type",x),1], age = spwagevals["age",1]*spwagereg_agemult,agesq=spwagevals["age2",1]/100, health=c(0,spwagevals["mod",1]))
                  )),
                  agelength = agelength
    )
    
  }
  if(post== FALSE & covs == FALSE){
    wagepars<-list( wagevar= (0.0251664^0.5)/(agelength^0.5),
                    wagecorr= 0.206,
                    wagereg = list(list(cons= -0.019831+shift*(0.0251664^0.5)/(agelength^0.5)+ list(0.7197562,1.279714,1.78743)[[1]],age = 0.0534633*wagereg_agemult,agesq=-0.000516491,health=c(0, -0.0230756,-0.1113683)),
                                   list(cons= -0.019831+shift*(0.0251664^0.5)/(agelength^0.5)+ list(0.7197562,1.279714,1.78743)[[2]],age = 0.0534633*wagereg_agemult,agesq=-0.000516491,health=c(0,-0.0230756,-0.1113683)),
                                   list(cons= -0.019831+shift*(0.0251664^0.5)/(agelength^0.5)+ list(0.7197562,1.279714,1.78743)[[3]],age = 0.0534633*wagereg_agemult,agesq=-0.000516491,health=c(0,-0.0230756,-0.1113683))
                    ),
                    agelength = agelength
    )
    
    spwagepars<-list( wagevar= (0.0582273^0.5)/(agelength^0.5),
                      wagecorr= 0.206,
                      wagereg = list(list(cons= -0.3255746 +spshift*(0.0582273^0.5)/(agelength^0.5)+list(1.25655 ,1.549101, 1.72255)[[1]],age = 0.0462337*spwagereg_agemult,agesq=-0.000427691,health=c(0,-0.0243041)),
                                     list(cons= -0.3255746 +spshift*(0.0582273^0.5)/(agelength^0.5)+list(1.25655 ,1.549101, 1.72255)[[2]],age = 0.0462337*spwagereg_agemult,agesq=-0.000427691,health=c(0,-0.0243041)),
                                     list(cons= -0.3255746 +spshift*(0.0582273^0.5)/(agelength^0.5)+list(1.25655 ,1.549101, 1.72255)[[3]],age = 0.0462337*spwagereg_agemult,agesq=-0.000427691,health=c(0,-0.0243041))
                      ),
                      agelength = agelength
    )
  }
  if(post == TRUE & covs == FALSE){
    wagepars<-list( wagevar= (0.0242^0.5)/(agelength^0.5),
                    wagecorr= 0.153,
                    wagereg = list(list(cons= shift*0.15939824/(agelength^0.5)+ list(0.504,1.166,1.761)[[1]],age = 0.0586*wagereg_agemult,agesq=-0.000563,health=c(0, -0.0199,-0.0778)),
                                   list(cons= shift*0.15939824/(agelength^0.5)+ list(0.504,1.166,1.761)[[2]],age = 0.0586*wagereg_agemult,agesq=-0.000563,health=c(0,-0.0199,-0.0778)),
                                   list(cons= shift*0.15939824/(agelength^0.5)+ list(0.504,1.166,1.761)[[3]],age = 0.0586*wagereg_agemult,agesq=-0.000563,health=c(0,-0.0199,-0.0778))
                    ),
                    agelength = agelength
    )
    
    spwagepars<-list( wagevar= (0.0396^0.5)/(agelength^0.5),
                      wagecorr= 0.153,
                      wagereg = list(list(cons= spshift*0.19688245/(agelength^0.5)+list(1.464 ,1.698, 1.918)[[1]],age = 0.0316*spwagereg_agemult,agesq=-0.000368,health=c(0,-0.00131)),
                                     list(cons= spshift*0.19688245/(agelength^0.5)+list(1.464 ,1.698, 1.918)[[2]],age = 0.0316*spwagereg_agemult,agesq=-0.000368,health=c(0,-0.00131)),
                                     list(cons= spshift*0.19688245/(agelength^0.5)+list(1.464 ,1.698, 1.918)[[3]],age = 0.0316*spwagereg_agemult,agesq=-0.000368,health=c(0,-0.00131))
                      ),
                      agelength = agelength
    )
    
  }
  Wage<-array(NA,c(numtypes,50*agelength,wage_gridsize))
  spWage<-array(NA,c(numtypes,50*agelength,spwage_gridsize))
  
  for(t in 1:50*agelength){
    gridpoints<-qlnorm(seq(0.01,0.99,length.out=wage_gridsize),meanlog=0,sdlog=wagepars$wagevar*(t)^0.5)
    if(spwage_gridsize>1) spgridpoints<-qlnorm(seq(0.05,0.95,length.out=spwage_gridsize),meanlog=0,sdlog=spwagepars$wagevar*(t)^0.5)
    else spgridpoints <- qlnorm(0.5,meanlog=0,sdlog=spwagepars$wagevar*(t)^0.5)
    for(type in 1:numtypes){
      Wage[type,t,]<-array(PermWage(t+23*wagepars$agelength,wagepars,type),dim=wage_gridsize)
      Wage[type,t,]<- Wage[type,t,]*gridpoints
      spWage[type,t,]<-array(PermWage(t+23*spwagepars$agelength,spwagepars,type),dim=spwage_gridsize)
      spWage[type,t,]<- spWage[type,t,]*spgridpoints
    }
  }
  return(list(wagepars=wagepars,spwagepars=spwagepars,Wage=Wage, spWage=spWage))
}



wagesetup <- function(shift=0,spshift=0,wage_gridsize=10, spwage_gridsize=10,
                      wagereg_agemult=1,
                      spwagereg_agemult=1,covs=TRUE,post=FALSE,filesuffix=""){
  wagevals<-read.table(paste0("wageparams",filesuffix,".txt"),row.names=1,fill=TRUE)
  spwagevals<-read.table(paste0("spwageparams",filesuffix,".txt"),row.names=1,fill=TRUE)
  numtypes<-sum(grepl("Type",rownames(wagevals)))

  wagevars<-read.table(paste0("wagevarcov",filesuffix,".txt"),row.names=1)
  
  if(post==FALSE & covs == TRUE){
    wagepars <- list(wagevar = wagevars[1,1]^0.5/(agelength^0.5),
                     wagecorr = wagevars[3,1]
                     )
    wagepars <-c(wagepars,
                 wagereg = list(lapply(1:numtypes,
                                  function(x) list(cons = 0+shift*wagepars$wagevar+wagevals[paste0("Type",x),1], age = wagevals["diffage",1]*wagereg_agemult,agesq=wagevals["diffage2",1]/100, health=c(0,wagevals["diffmod",1],wagevals["diffsev",1]))
                                  )),
                 agelength = agelength
              )
    
    spwagepars <- list(wagevar = wagevars[2,1]^0.5/(agelength^0.5),
                       wagecorr = wagevars[3,1]
    )
    spwagepars<-c(spwagepars,
                  wagereg = list(lapply(1:numtypes,
                                        function(x) list(cons = 0+spshift*spwagepars$wagevar+spwagevals[paste0("Type",x),1], age = spwagevals["diffage",1]*spwagereg_agemult,agesq=spwagevals["diffage2",1]/100, health=c(0,spwagevals["diffmod",1]))
                  )),
                  agelength = agelength
                  )

  }
  if(post== FALSE & covs == FALSE){
    wagepars<-list( wagevar= (0.0251664^0.5)/(agelength^0.5),
                    wagecorr= 0.206,
                    wagereg = list(list(cons= -0.019831+shift*(0.0251664^0.5)/(agelength^0.5)+ list(0.7197562,1.279714,1.78743)[[1]],age = 0.0534633*wagereg_agemult,agesq=-0.000516491,health=c(0, -0.0230756,-0.1113683)),
                                   list(cons= -0.019831+shift*(0.0251664^0.5)/(agelength^0.5)+ list(0.7197562,1.279714,1.78743)[[2]],age = 0.0534633*wagereg_agemult,agesq=-0.000516491,health=c(0,-0.0230756,-0.1113683)),
                                   list(cons= -0.019831+shift*(0.0251664^0.5)/(agelength^0.5)+ list(0.7197562,1.279714,1.78743)[[3]],age = 0.0534633*wagereg_agemult,agesq=-0.000516491,health=c(0,-0.0230756,-0.1113683))
                    ),
                    agelength = agelength
    )
    
    spwagepars<-list( wagevar= (0.0582273^0.5)/(agelength^0.5),
                      wagecorr= 0.206,
                      wagereg = list(list(cons= -0.3255746 +spshift*(0.0582273^0.5)/(agelength^0.5)+list(1.25655 ,1.549101, 1.72255)[[1]],age = 0.0462337*spwagereg_agemult,agesq=-0.000427691,health=c(0,-0.0243041)),
                                     list(cons= -0.3255746 +spshift*(0.0582273^0.5)/(agelength^0.5)+list(1.25655 ,1.549101, 1.72255)[[2]],age = 0.0462337*spwagereg_agemult,agesq=-0.000427691,health=c(0,-0.0243041)),
                                     list(cons= -0.3255746 +spshift*(0.0582273^0.5)/(agelength^0.5)+list(1.25655 ,1.549101, 1.72255)[[3]],age = 0.0462337*spwagereg_agemult,agesq=-0.000427691,health=c(0,-0.0243041))
                      ),
                      agelength = agelength
    )
  }
  if(post == TRUE & covs == FALSE){
    wagepars<-list( wagevar= (0.0242^0.5)/(agelength^0.5),
                    wagecorr= 0.153,
                    wagereg = list(list(cons= shift*0.15939824/(agelength^0.5)+ list(0.504,1.166,1.761)[[1]],age = 0.0586*wagereg_agemult,agesq=-0.000563,health=c(0, -0.0199,-0.0778)),
                                   list(cons= shift*0.15939824/(agelength^0.5)+ list(0.504,1.166,1.761)[[2]],age = 0.0586*wagereg_agemult,agesq=-0.000563,health=c(0,-0.0199,-0.0778)),
                                   list(cons= shift*0.15939824/(agelength^0.5)+ list(0.504,1.166,1.761)[[3]],age = 0.0586*wagereg_agemult,agesq=-0.000563,health=c(0,-0.0199,-0.0778))
                    ),
                    agelength = agelength
    )
    
    spwagepars<-list( wagevar= (0.0396^0.5)/(agelength^0.5),
                      wagecorr= 0.153,
                      wagereg = list(list(cons= spshift*0.19688245/(agelength^0.5)+list(1.464 ,1.698, 1.918)[[1]],age = 0.0316*spwagereg_agemult,agesq=-0.000368,health=c(0,-0.00131)),
                                     list(cons= spshift*0.19688245/(agelength^0.5)+list(1.464 ,1.698, 1.918)[[2]],age = 0.0316*spwagereg_agemult,agesq=-0.000368,health=c(0,-0.00131)),
                                     list(cons= spshift*0.19688245/(agelength^0.5)+list(1.464 ,1.698, 1.918)[[3]],age = 0.0316*spwagereg_agemult,agesq=-0.000368,health=c(0,-0.00131))
                      ),
                      agelength = agelength
    )
    
  }
  Wage<-array(NA,c(numtypes,50*agelength,wage_gridsize))
  spWage<-array(NA,c(numtypes,50*agelength,spwage_gridsize))
  
  for(t in 1:50*agelength){
    gridpoints<-qlnorm(seq(0.01,0.99,length.out=wage_gridsize),meanlog=0,sdlog=wagepars$wagevar*(t)^0.5)
    if(spwage_gridsize>1) spgridpoints<-qlnorm(seq(0.05,0.95,length.out=spwage_gridsize),meanlog=0,sdlog=spwagepars$wagevar*(t)^0.5)
    else spgridpoints <- qlnorm(0.5,meanlog=0,sdlog=spwagepars$wagevar*(t)^0.5)
    for(type in 1:numtypes){
      Wage[type,t,]<-array(PermWage(t+23*wagepars$agelength,wagepars,type),dim=wage_gridsize)
      Wage[type,t,]<- Wage[type,t,]*gridpoints
      spWage[type,t,]<-array(PermWage(t+23*spwagepars$agelength,spwagepars,type),dim=spwage_gridsize)
      spWage[type,t,]<- spWage[type,t,]*spgridpoints
    }
  }
  return(list(wagepars=wagepars,spwagepars=spwagepars,Wage=Wage, spWage=spWage))
}


wagesetup_jmp <- function(shift=0,spshift=0,wage_gridsize=10, spwage_gridsize=10,
                      wagereg_agemult=1,
                      spwagereg_agemult=1,post=FALSE){
  if(post==FALSE){
  wagepars<-list( wagevar= (0.0234761^0.5)/(agelength^0.5),
                  wagecorr= 0.18367903,
                  wagereg = list(list(cons= -0.0235764+shift*0.15939824/(agelength^0.5)+ list(0.6094251,1.159559,1.668392)[[1]],age = 0.0587180*wagereg_agemult,agesq=-0.00056567,health=c(0, -0.0207469,-0.0786841)),
                                 list(cons= -0.0235764+shift*0.15939824/(agelength^0.5)+ list(0.6094251,1.159559,1.668392)[[2]],age = 0.0587180*wagereg_agemult,agesq=-0.00056567,health=c(0,-0.0207469,-0.0786841)),
                                 list(cons= -0.0235764+shift*0.15939824/(agelength^0.5)+ list(0.6094251,1.159559,1.668392)[[3]],age = 0.0587180*wagereg_agemult,agesq=-0.00056567,health=c(0,-0.0207469,-0.0786841))
                  ),
                  agelength = agelength
  )
  
  spwagepars<-list( wagevar= (0.0399729^0.5)/(agelength^0.5),
                    wagecorr= 0.18367903,
                    wagereg = list(list(cons= 0.1973193 +spshift*0.19688245/(agelength^0.5)+list(1.337045 ,1.535034, 1.733289)[[1]],age = 0.0300251*spwagereg_agemult,agesq=-0.000362589,health=c(0,-0.0055049)),
                                   list(cons= 0.1973193 +spshift*0.19688245/(agelength^0.5)+list(1.337045 ,1.535034, 1.733289)[[2]],age = 0.0300251*spwagereg_agemult,agesq=-0.000362589,health=c(0,-0.0055049)),
                                   list(cons= 0.1973193 +spshift*0.19688245/(agelength^0.5)+list(1.337045 ,1.535034, 1.733289)[[3]],age = 0.0300251*spwagereg_agemult,agesq=-0.000362589,health=c(0,-0.0055049))
                    ),
                    agelength = agelength
  )
  }
  else{
    wagepars<-list( wagevar= (0.0242^0.5)/(agelength^0.5),
                    wagecorr= 0.153,
                    wagereg = list(list(cons= shift*0.15939824/(agelength^0.5)+ list(0.504,1.166,1.761)[[1]],age = 0.0586*wagereg_agemult,agesq=-0.000563,health=c(0, -0.0199,-0.0778)),
                                   list(cons= shift*0.15939824/(agelength^0.5)+ list(0.504,1.166,1.761)[[2]],age = 0.0586*wagereg_agemult,agesq=-0.000563,health=c(0,-0.0199,-0.0778)),
                                   list(cons= shift*0.15939824/(agelength^0.5)+ list(0.504,1.166,1.761)[[3]],age = 0.0586*wagereg_agemult,agesq=-0.000563,health=c(0,-0.0199,-0.0778))
                    ),
                    agelength = agelength
    )
    
    spwagepars<-list( wagevar= (0.0396^0.5)/(agelength^0.5),
                      wagecorr= 0.153,
                      wagereg = list(list(cons= spshift*0.19688245/(agelength^0.5)+list(1.464 ,1.698, 1.918)[[1]],age = 0.0316*spwagereg_agemult,agesq=-0.000368,health=c(0,-0.00131)),
                                     list(cons= spshift*0.19688245/(agelength^0.5)+list(1.464 ,1.698, 1.918)[[2]],age = 0.0316*spwagereg_agemult,agesq=-0.000368,health=c(0,-0.00131)),
                                     list(cons= spshift*0.19688245/(agelength^0.5)+list(1.464 ,1.698, 1.918)[[3]],age = 0.0316*spwagereg_agemult,agesq=-0.000368,health=c(0,-0.00131))
                      ),
                      agelength = agelength
    )
    
  }
  Wage<-array(NA,c(numtypes,50*agelength,wage_gridsize))
  spWage<-array(NA,c(numtypes,50*agelength,spwage_gridsize))
  
  for(t in 1:50*agelength){
    gridpoints<-qlnorm(seq(0.01,0.99,length.out=wage_gridsize),meanlog=0,sdlog=wagepars$wagevar*(t)^0.5)
    if(spwage_gridsize>1) spgridpoints<-qlnorm(seq(0.01,0.99,length.out=spwage_gridsize),meanlog=0,sdlog=spwagepars$wagevar*(t)^0.5)
    else spgridpoints <- qlnorm(0.5,meanlog=0,sdlog=spwagepars$wagevar*(t)^0.5)
    for(type in 1:numtypes){
      Wage[type,t,]<-array(PermWage(t+23*wagepars$agelength,wagepars,type),dim=wage_gridsize)
      Wage[type,t,]<- Wage[type,t,]*gridpoints
      spWage[type,t,]<-array(PermWage(t+23*spwagepars$agelength,spwagepars,type),dim=spwage_gridsize)
      spWage[type,t,]<- spWage[type,t,]*spgridpoints
    }
  }
  return(list(wagepars=wagepars,spwagepars=spwagepars,Wage=Wage, spWage=spWage))
}

#original LP equation stuff:
wagesetup_lp <- function(shift=0,spshift=0,wage_gridsize=10){
  wagepars<-list( wagevar= 0.1632008/(agelength^0.5),
                  wagereg = list(list(cons= -0.0314009+shift*0.1632008/(agelength^0.5)+ list(1.03208,1.593764,2.185418)[[1]],age = 0.0524131,agesq=-0.000671353,health=c(0, -0.0568435,-0.1771249)),
                                 list(cons= -0.0314009+shift*0.1632008/(agelength^0.5)+ list(1.03208,1.593764,2.185418)[[2]],age = 0.0524131,agesq=-0.000671353,health=c(0,-0.0568435,-0.1771249)),
                                 list(cons= -0.0314009+shift*0.1632008/(agelength^0.5)+ list(1.03208,1.593764,2.185418)[[3]],age = 0.0524131,agesq=-0.000671353,health=c(0,-0.0568435,-0.1771249))
                  ),
                  agelength = agelength
  )
  
  spwagepars<-list( wagevar= 0.182/(agelength^0.5),
                    wagereg = list(list(cons= -0.2134809+spshift*0.182/(agelength^0.5)+list(0.656225 ,0.7970078,0.9082671)[[1]],age = 0.0620146,agesq=-0.0006555559,health=c(0,0.0037686)),
                                   list(cons= -0.2134809+spshift*0.182/(agelength^0.5)+list(0.656225 ,0.7970078,0.9082671)[[2]],age = 0.0620146,agesq=-0.0006555559,health=c(0,0.0037686)),
                                   list(cons= -0.2134809+spshift*0.182/(agelength^0.5)+list(0.656225 ,0.7970078,0.9082671)[[3]],age = 0.0620146,agesq=-0.0006555559,health=c(0,0.0037686))
                    ),
                    agelength = agelength
  )
  ##this is from wage regression, no restriction on birthy and no imputation of missing spousal fixed effects.
  # spwagepars<-list( wagevar= 0.182/(agelength^0.5),
  #                   wagereg = list(list(cons= -0.2134809+spshift*0.182/(agelength^0.5)+list(0.8742635 ,1.108983,1.284846)[[1]],age = 0.0704852,agesq=-0.000670589,health=c(0,-0.0502703)),
  #                                  list(cons= -0.2134809+spshift*0.182/(agelength^0.5)+list(0.8742635 ,1.108983,1.284846)[[2]],age = 0.0704852,agesq=-0.000670589,health=c(0,-0.0502703)),
  #                                  list(cons= -0.2134809+spshift*0.182/(agelength^0.5)+list(0.8742635 ,1.108983,1.284846)[[3]],age = 0.0704852,agesq=-0.000670589,health=c(0,-0.0502703))
  #                   ))  
  # #spousal wage eq, restricting head to >=1945, assigning missing values to 10th %tile
  # spwagepars<-list( wagevar= 0.182/(agelength^0.5),
  #                   wagereg = list(list(cons= 0.0317059+spshift*0.182/(agelength^0.5)+list(0.3278311 ,0.4480985,0.5557425)[[1]],age = 0.0696291,agesq=-0.000771306,health=c(0,0.0172439)),
  #                                  list(cons= 0.0317059+spshift*0.182/(agelength^0.5)+list(0.3278311 ,0.4480985,0.5557425)[[2]],age = 0.0696291,agesq=-0.000771306,health=c(0,0.0172439)),
  #                                  list(cons= 0.0317059+spshift*0.182/(agelength^0.5)+list(0.3278311 ,0.4480985,0.5557425)[[3]],age = 0.0696291,agesq=-0.000771306,health=c(0,0.0172439))
  #                   ),
  #                   agelength = agelength
  # )
  
  Wage<-array(NA,c(numtypes,50*agelength,wage_gridsize))
  
  for(t in 1:50*agelength){
    gridpoints<-qlnorm(seq(0.01,0.99,length.out=wage_gridsize),meanlog=0,sdlog=wagepars$wagevar*(t)^0.5)
    
    for(type in 1:numtypes){
      Wage[type,t,]<-array(PermWage(t+23*wagepars$agelength,wagepars,type),dim=wage_gridsize)
      Wage[type,t,]<- Wage[type,t,]*gridpoints
    }
  }
  return(list(wagepars=wagepars,spwagepars=spwagepars,Wage=Wage))
}

#LP stuff, my data:
wagesetup_lp2 <- function(shift=0,spshift=0,wage_gridsize=10){
  wagepars<-list( wagevar= 0.1632008/(agelength^0.5),
                  wagereg = list(list(cons= -0.0362812+shift*0.1632008/(agelength^0.5)+ list(0.6268305,1.219131,1.772136)[[1]],age = 0.0562602,agesq=-0.000551751,health=c(0, -0.0191509,-0.0669418)),
                                 list(cons= -0.0362812+shift*0.1632008/(agelength^0.5)+ list(0.6268305,1.219131,1.772136)[[2]],age = 0.0562602,agesq=-0.000551751,health=c(0,-0.0191509,-0.0669418)),
                                 list(cons= -0.0362812+shift*0.1632008/(agelength^0.5)+ list(0.6268305,1.219131,1.772136)[[3]],age = 0.0562602,agesq=-0.000551751,health=c(0,-0.0191509,-0.0669418))
                  ),
                  agelength = agelength
  )
  
  spwagepars<-list( wagevar= 0.182/(agelength^0.5),
                    wagereg = list(list(cons= 0.1224924+spshift*0.182/(agelength^0.5)+list(0.8879625 ,1.045816,1.168454)[[1]],age = 0.04,agesq=-0.000425044,health=c(0,0.00312)),
                                   list(cons= 0.1224924+spshift*0.182/(agelength^0.5)+list(0.8879625 ,1.045816,1.168454)[[2]],age = 0.04,agesq=-0.000425044,health=c(0,0.00312)),
                                   list(cons= 0.1224924+spshift*0.182/(agelength^0.5)+list(0.8879625 ,1.045816,1.168454)[[3]],age = 0.04,agesq=-0.000425044,health=c(0,0.00312))
                    ),
                    agelength = agelength
  )
  
  
  Wage<-array(NA,c(numtypes,50*agelength,wage_gridsize))
  
  for(t in 1:50*agelength){
    gridpoints<-qlnorm(seq(0.01,0.99,length.out=wage_gridsize),meanlog=0,sdlog=wagepars$wagevar*(t)^0.5)
    
    for(type in 1:numtypes){
      Wage[type,t,]<-array(PermWage(t+23*wagepars$agelength,wagepars,type),dim=wage_gridsize)
      Wage[type,t,]<- Wage[type,t,]*gridpoints
    }
  }
  return(list(wagepars=wagepars,spwagepars=spwagepars,Wage=Wage))
}

statesetup<-function(asset_gridsize=20,
                     wage_gridsize=10, marginal=1,
                    marriagehealth=1,dir,filesuffix=""){
  origdir<-getwd()
  setwd(dir)
  if(marginal==1) temp<-as.data.table(read.csv(paste0("empiricalhealthprobs_",filesuffix,"_marginal.csv"),row.names=NULL))
  else temp<-as.data.table(read.csv(paste0("empiricalhealthprobs_",filesuffix,".csv"),row.names=NULL))
  tempfam_mar<-as.data.table(read.csv(paste0("famsize_married",filesuffix,".csv"),row.names=NULL))
  tempfam_sing<-as.data.table(read.csv(paste0("famsize_single",filesuffix,".csv"),row.names=NULL))
  
  if(marriagehealth==1){
  tempmar<-as.data.table(read.csv(paste0("empiricalmarriageprobs_healthy",filesuffix,".csv"),row.names=NULL))
  tempmarstock<-as.data.table(read.csv(paste0("empiricalmarriagestocks_healthy",filesuffix,".csv"),row.names=NULL))
  }
  else{
    tempmar<-as.data.table(read.csv(paste0("empiricalmarriageprobs_",filesuffix,".csv"),row.names=NULL))
    tempmarstock<-as.data.table(read.csv(paste0("empiricalmarriagestocks_",filesuffix,".csv"),row.names=NULL))
  }
  
  assetdist<-as.data.table(read.csv("assetdist_byage.csv",row.names=NULL))
  
  setwd(origdir)
  
  healthprob<-list()
  if(marginal==1){
    for(typ in sort(unique(temp$type))){
      healthprob[[typ]]<-array(NA,dim=c(50*agelength,6,6))
      for(s in 1:3){
        for(sp in 1:2){
          pos<-1
          for(t in 1:50){
            for(rep in 1:agelength){
              healthprob[[typ]][pos,(s-1)*2+sp,]<-as.vector(t(outer(temp[type==typ&from==s & age == t+22,p]
                                                                    ,temp[type==typ&from==sp & age == t+22 & to !=3,sp_p]))) #probabilities  sum to 1
              pos<-pos+1
            }
          }
        }
      }
    }
  }
  else{
    for(typ in sort(unique(temp$type))){
      healthprob[[typ]]<-array(NA,dim=c(50*agelength,6,6))
      for(s in 1:3){
        for(sp in 1:2){
          pos<-1
          for(t in 1:50){
            for(rep in 1:agelength){
              healthprob[[typ]][pos,(s-1)*2+sp,]<-temp[type==typ&from==(s-1)*2+sp & age == t+22,p] #L&P normalize by making them sum to 1
              pos<-pos+1
            }
          }
        }
      }
    }
  }
  #adjusting transition probabilities:(1/agelength) prob. of transitioning
  for(typ in sort(unique(temp$type))){
    for(t in 1:50*agelength){
      for(s in 1:6){
        healthprob[[typ]][t,s,] = healthprob[[typ]][t,s,]*(1/agelength)
        healthprob[[typ]][t,s,s] = healthprob[[typ]][t,s,s]+(1-1/agelength)
      }
    }
  }
  if(marriagehealth==1){
  marriageprob<-list()
  for(typ in sort(unique(temp$type))){
    marriageprob[[typ]]<-list()
    marriageprob[[typ]][[1]]<-array(NA,dim=c(50*agelength,2))
    marriageprob[[typ]][[2]]<-array(NA,dim=c(50*agelength,2))
    marriageprob[[typ]][[3]]<-array(NA,dim=c(50*agelength,2))
    for(m in 1:2){
      pos<-1
      for(t in 1:50){
        for(rep in 1:agelength){
          marriageprob[[typ]][[1]][pos,3-m]<-tempmar[type==typ&from==m & age == t+22 & to ==2,p] #L&P normalize by making them sum to 1
          pos<-pos+1
        }
      }
    }
    marriageprob[[typ]][[2]]<-marriageprob[[typ]][[1]]
    marriageprob[[typ]][[3]]<-marriageprob[[typ]][[1]]
    marriageprob[[typ]][[3]][,2]<-pmax(marriageprob[[typ]][[1]][,2] - 0,0) #STOPPED THIS ADJUSTMENT #severely disabled singles are less likely to get married, per event studies
  }
  #adjusting transition probabilities:(1/agelength) prob. of transitioning
  for(h in 1:3){
  for(typ in sort(unique(temp$type))){
    for(t in 1:50*agelength){
      marriageprob[[typ]][[h]][t,] = marriageprob[[typ]][[h]][t,]*(1/agelength)
      marriageprob[[typ]][[h]][t,1] = marriageprob[[typ]][[h]][t,1]+(1-1/agelength)
    }
  }
  }
  }
  else{
    marriageprob<-list()
    for(typ in sort(unique(temp$type))){
      marriageprob[[typ]]<-array(NA,dim=c(50*agelength,2))
      for(m in 1:2){
        pos<-1
        for(t in 1:50){
          for(rep in 1:agelength){
            marriageprob[[typ]][pos,3-m]<-tempmar[type==typ&from==m & age == t+22 & to ==2,p] #L&P normalize by making them sum to 1
            pos<-pos+1
          }
        }
      }
    }
    #adjusting transition probabilities:(1/agelength) prob. of transitioning
    for(typ in sort(unique(temp$type))){
      for(t in 1:50*agelength){
        marriageprob[[typ]][t,] = marriageprob[[typ]][t,]*(1/agelength)
        marriageprob[[typ]][t,1] = marriageprob[[typ]][t,1]+(1-1/agelength)
      }
    }
  }
  
  marriageprob_init <- tempmarstock[age==23,M]
  
  numadults<-list()
  numkids<-list()
  for(typ in sort(unique(temp$type))){
    numadults[[typ]]<-array(NA,dim=c(50*agelength,2))
    numkids[[typ]]<-array(NA,dim=c(50*agelength,2))
  }
  
  pos<-1
  for(t in 1:50){
    for(rep in 1:agelength){
      for(typ in sort(unique(temp$type))){
        if(t < 50){
          numadults[[typ]][pos,1]<-tempfam_mar[type_prodwage==typ & age==t+22,numadults_p]*(agelength-rep+1)/agelength+tempfam_mar[type_prodwage==typ&age==t+23,numadults_p]*(rep-1)/agelength
          numadults[[typ]][pos,2]<-tempfam_sing[type_prodwage==typ & age==t+22,numadults_p]*(agelength-rep+1)/agelength+tempfam_sing[type_prodwage==typ&age==t+23,numadults_p]*(rep-1)/agelength
          
          numkids[[typ]][pos,1]<-tempfam_mar[type_prodwage==typ&age==t+22,kids_p]*(agelength-rep+1)/agelength+tempfam_mar[type_prodwage==typ&age==t+23,kids_p]*(rep-1)/agelength
          numkids[[typ]][pos,2]<-tempfam_sing[type_prodwage==typ&age==t+22,kids_p]*(agelength-rep+1)/agelength+tempfam_sing[type_prodwage==typ&age==t+23,kids_p]*(rep-1)/agelength
        } else {
          numadults[[typ]][pos,1]<-tempfam_mar[(t-1)*length(unique(temp$type))+typ,numadults_p]
        numadults[[typ]][pos,2]<-tempfam_sing[(t-1)*length(unique(temp$type))+typ,numadults_p]
        numkids[[typ]][pos,1]<-tempfam_mar[(t-1)*length(unique(temp$type))+typ,kids_p]
        numkids[[typ]][pos,2]<-tempfam_sing[(t-1)*length(unique(temp$type))+typ,kids_p]
        }
      }
      pos<-pos+1
    }
  }
  numkids<-lapply(numkids,function(x) pmax(x,0))
  numadults<-lapply(numadults,function(x) pmax(x,0))
  foodstamps<-list()
  for(typ in sort(unique(temp$type))){
    foodstamps[[typ]]<-array(NA,dim=c(50*agelength,2))
  }
  
  for(typ in sort(unique(temp$type))){
    tot<-numadults[[typ]]+numkids[[typ]]
    for(t in 1:nrow(tot)){
      for(m in 1:2){
        if(tot[t,m] <=1) foodstamps[[typ]][t,m]<-119*tot[t,m]*12/agelength
        if(tot[t,m] > 1 & tot[t,m]<=2) foodstamps[[typ]][t,m]<- (119+(218-119)*(tot[t,m]-1))*12/agelength
        if(tot[t,m] > 2 & tot[t,m]<=3) foodstamps[[typ]][t,m] <- (218+(313-218)*(tot[t,m]-2))*12/agelength
        if(tot[t,m] > 3) foodstamps[[typ]][t,m] <- (313+(397-313)*(tot[t,m]-3))*12/agelength
      }
    }
  }
  
  
  oecd<-list()
  for(x in sort(unique(temp$type))){
    oecd[[x]]<-1+(numadults[[x]]-1)*0.5 + numkids[[x]]*0.3
  }
  
  #EMPIRICALLY DRIVEN ASSET GRID:
  Asset<-array(NA,c(numtypes,50*agelength+1,asset_gridsize))
  for (typ in sort(unique(temp$type))){
    assetlow<-rep(0,50*agelength+1)
    for(t in 1:40*agelength){
      Asset[typ,t,]<-c(assetlow[t],assetdist[type_fix==typ&age==t+22&p%in%((1:(asset_gridsize-1))*100/asset_gridsize),asset])
    }
    for(t in (40*agelength+1):(50*agelength+1)){
      Asset[typ,t,]<-Asset[typ,t-1,]
    }
    #MAKING SURE GRID IS MONOTONICALLY INCREASING:
    for(t in 1:(50*agelength+1)){
      for(p in (asset_gridsize-1):2){
        if(Asset[typ,t,p] <= Asset[typ,t,p-1] | Asset[typ,t,p]>=Asset[typ,t,p+1]) Asset[typ,t,p]<-0.5*Asset[typ,t,p-1]+0.5*Asset[typ,t,p+1]
      }
    }
  }
  #health and marriage probabilities need to be shifted ahead 1 period:
  healthprob<-lapply(healthprob,function(x) x[-1,,])
  if(marriagehealth==0) marriageprob<-lapply(marriageprob,function(x) x[-1,])
  else marriageprob<-lapply(marriageprob,function(x) lapply(x, function(z) z[-1,]))
  states<-list(Asset=Asset,foodstamps=foodstamps,numkids=numkids,healthprob=healthprob,oecd=oecd, marriageprob=marriageprob,marriageprob_init=marriageprob_init)
  return(states)
}

statesetup_old<-function(asset_gridsize=20,
                         wage_gridsize=10,dir){
  origdir<-getwd()
  setwd(dir)
  temp<-as.data.table(read.csv("empiricalhealthprobs_married.csv",row.names=NULL))
  tempfam<-as.data.table(read.csv("famsize_married.csv",row.names=NULL))
  
  setwd(origdir)
  
  assetdist<-as.data.table(read.csv("assetdist_byage.csv",row.names=NULL))
  
  
  healthprob<-list()
  for(typ in 1:numtypes){
    healthprob[[typ]]<-array(NA,dim=c(50*agelength,6,6))
    for(s in 1:3){
      for(sp in 1:2){
        pos<-1
        for(t in 1:50){
          for(rep in 1:agelength){
            healthprob[[typ]][pos,(s-1)*2+sp,]<-temp[type==typ&from==(s-1)*2+sp & age == t+22,p] #L&P normalize by making them sum to 1
            pos<-pos+1
          }
        }
      }
    }
  }
  #adjusting transition probabilities:(1/agelength) prob. of transitioning
  for(typ in 1:numtypes){
    for(t in 1:50*agelength){
      for(s in 1:6){
        healthprob[[typ]][t,s,] = healthprob[[typ]][t,s,]*(1/agelength)
        healthprob[[typ]][t,s,s] = healthprob[[typ]][t,s,s]+(1-1/agelength)
      }
    }
  }
  
  
  
  
  numadults<-list()
  numkids<-list()
  for(typ in 1:numtypes){
    numadults[[typ]]<-array(NA,dim=c(50*agelength,1))
    numkids[[typ]]<-array(NA,dim=c(50*agelength,1))
  }
  
  pos<-1
  for(t in 1:50){
    for(rep in 1:agelength){
      for(typ in 1:numtypes){
        if(t < 50){
          numadults[[typ]][pos,]<-tempfam[type_fix==typ & age==t+22,numadults_p]*(agelength-rep+1)/agelength+tempfam[type_fix==typ&age==t+23,numadults_p]*(rep-1)/agelength
          numkids[[typ]][pos,]<-tempfam[type_fix==typ&age==t+22,kids_p]*(agelength-rep+1)/agelength+tempfam[type_fix==typ&age==t+23,kids_p]*(rep-1)/agelength
        }
        else
          numadults[[typ]][pos,]<-tempfam[(t-1)*numtypes+typ,numadults_p]
        numkids[[typ]][pos,]<-tempfam[(t-1)*numtypes+typ,kids_p]
        
      }
      pos<-pos+1
    }
  }
  foodstamps<-list()
  for(typ in 1:numtypes){
    foodstamps[[typ]]<-array(NA,dim=c(50*agelength,1))
  }
  
  for(typ in 1:numtypes){
    tot<-numadults[[typ]]+numkids[[typ]]
    for(t in 1:length(tot)){
      if(tot[t] <=1) foodstamps[[typ]][t]<-119*tot[t]*12/agelength
      if(tot[t] > 1 & tot[t]<=2) foodstamps[[typ]][t]<- (119+(218-119)*(tot[t]-1))*12/agelength
      if(tot[t] > 2 & tot[t]<=3) foodstamps[[typ]][t] <- (218+(313-218)*(tot[t]-2))*12/agelength
      if(tot[t] > 3) foodstamps[[typ]][t] <- (313+(397-313)*(tot[t]-3))*12/agelength
    }
  }
  
  
  oecd<-list()
  for(x in 1:numtypes){
    oecd[[x]]<-1+(numadults[[x]]-1)*0.5 + numkids[[x]]*0.3
  }
  
  #EMPIRICALLY DRIVEN ASSET GRID:
  Asset<-array(NA,c(numtypes,50*agelength+1,asset_gridsize))
  for (typ in 1:numtypes){
    assetlow<-rep(0,50*agelength+1)
    for(t in 1:40*agelength){
      Asset[typ,t,]<-c(assetlow[t],assetdist[type_fix==typ&age==t+22&p%in%((1:(asset_gridsize-1))*100/asset_gridsize),asset])
    }
    for(t in (40*agelength+1):(50*agelength+1)){
      Asset[typ,t,]<-Asset[typ,t-1,]
    }
    #MAKING SURE GRID IS MONOTONICALLY INCREASING:
    for(t in 1:(50*agelength+1)){
      for(p in (asset_gridsize-1):2){
        if(Asset[typ,t,p] <= Asset[typ,t,p-1] | Asset[typ,t,p]>=Asset[typ,t,p+1]) Asset[typ,t,p]<-0.5*Asset[typ,t,p-1]+0.5*Asset[typ,t,p+1]
      }
    }
  }
  states<-list(Asset=Asset,foodstamps=foodstamps,numkids=numkids,healthprob=healthprob,oecd=oecd)
  return(states)
}


paramsetup<-function(paramvals,wagepars,spwagepars,
                     divorceshare=1,
                     health_halfdis=FALSE,
                     ptwork_halfcost=FALSE,
                     single_samecost=FALSE,
                     ptwork_halfdis=FALSE,
                     single_samedis=FALSE,
                     spouse_samejobrisk=FALSE,
                     single_samejobrisk=FALSE,
                     single_sameprefs = TRUE,
                     health_simpleprefs = TRUE,
                     single_sameallowanceprob=TRUE,
                     work_noallowanceprob=FALSE,
                     SSgenerosity=1,FSgenerosity=1,wagetaxmult=1,ficataxmult=1,
                     wageshifter=rep(1,50),spwageshifter=rep(1,50),
                     wagemult_range=c(1,1),
                     spwagemult_range=c(1,1),
                     goldsearch_tolerance=1e-3){
  
  params<-list(
    agelength=agelength,
    wagevar= matrix(c(wagepars$wagevar^2,wagepars$wagecorr*(wagepars$wagevar*spwagepars$wagevar),
                      wagepars$wagecorr*(wagepars$wagevar*spwagepars$wagevar),spwagepars$wagevar^2),2,2),
    #wagepars$wagevar,  #wage STDEV -- 0.027 #wagevar is variance of wage period shocks
    spwagevar = spwagepars$wagevar
    , tau = 1 #tau is the proportional tax used to fund social programs
    #THE FOLLOWING ARE PRE-SPECIFIED
    #    , reassess =  c(1-(1-0.222)^(4/agelength),1-(1-0.083)^(4/agelength),1-(1-0.036)^(4/agelength)) # reassess is the DI reassessment rate, vector,1 param per health level
    , reassess =  c(1,1/3,1/7) # reassess is the DI reassessment rate, vector,1 param per health level
    , divorceshare=divorceshare
    , gamma = 1.5
    , beta = 0.9756^(1/agelength)
    , R = 1.016^(1/agelength)
    ,wagereg_cons = c(0,0,0)
    ,wagereg_fe = c(wagepars$wagereg[[1]]$cons,wagepars$wagereg[[2]]$cons,wagepars$wagereg[[3]]$cons)
    ,wagereg_age1 = c(wagepars$wagereg[[1]]$age,wagepars$wagereg[[2]]$age,wagepars$wagereg[[3]]$age)
    ,wagereg_age2 = c(wagepars$wagereg[[1]]$agesq,wagepars$wagereg[[2]]$agesq,wagepars$wagereg[[3]]$agesq)
    ,wagereg_h1 = c(wagepars$wagereg[[1]]$health[1],wagepars$wagereg[[2]]$health[1],wagepars$wagereg[[2]]$health[1])
    ,wagereg_h2 = c(wagepars$wagereg[[1]]$health[2],wagepars$wagereg[[2]]$health[2],wagepars$wagereg[[2]]$health[2])
    ,wagereg_h3 = c(wagepars$wagereg[[1]]$health[3],wagepars$wagereg[[2]]$health[3],wagepars$wagereg[[2]]$health[3])
    ,spwagereg_cons = c(0,0,0)
    ,spwagereg_fe = c(spwagepars$wagereg[[1]]$cons,spwagepars$wagereg[[2]]$cons,spwagepars$wagereg[[3]]$cons)
    ,spwagereg_age1 = c(spwagepars$wagereg[[1]]$age,spwagepars$wagereg[[2]]$age,spwagepars$wagereg[[3]]$age)
    ,spwagereg_age2 = c(spwagepars$wagereg[[1]]$agesq,spwagepars$wagereg[[2]]$agesq,spwagepars$wagereg[[3]]$agesq)
    ,spwagereg_h1 = c(spwagepars$wagereg[[1]]$health[1],spwagepars$wagereg[[2]]$health[1],spwagepars$wagereg[[2]]$health[1])
    ,spwagereg_h2 = c(spwagepars$wagereg[[1]]$health[2],spwagepars$wagereg[[2]]$health[2],spwagepars$wagereg[[2]]$health[2]),
    hours=2000*agelength,
    kink1=5244/agelength, #formula cutoffs from 1996
    kink2=31620/agelength,
    kink3=62700/agelength,
    fsgross = 23920/agelength,
    fsnet= 18400/agelength, 
    UIcap = 12712/agelength,
    SSIval=5100/agelength,
    cost_reassess = 364,
    SSgenerosity=SSgenerosity,
    FSgenerosity=FSgenerosity,
    wagetaxmult=wagetaxmult,
    ficataxmult=ficataxmult,
    wageshifter=wageshifter,
    spwageshifter=spwageshifter,
    wagemults = runif(n=sum(numsims),min=min(wagemult_range),max=max(wagemult_range)),
    spwagemults = runif(n=sum(numsims),min=min(spwagemult_range),max=max(spwagemult_range)),
    goldsearch_tolerance=goldsearch_tolerance
  )
  
  F<-NULL
  F[1]<-paramvals["F1"]*params$hours*PermWage(23*agelength,wagepars,2)
  F[2]<-paramvals["F2"]*params$hours*PermWage(23*agelength,wagepars,2)
  F[3]<-paramvals["F3"]*params$hours*PermWage(23*agelength,wagepars,2)
  if(ptwork_halfcost==FALSE){
    F[4]<-paramvals["Fpart1"]*params$hours*PermWage(23*agelength,wagepars,2)
    F[5]<-paramvals["Fpart2"]*params$hours*PermWage(23*agelength,wagepars,2)
    F[6]<-paramvals["Fpart3"]*params$hours*PermWage(23*agelength,wagepars,2)
  }
  else{
    F[4]<-0/2
    F[5]<-paramvals["F2"]*params$hours*PermWage(23*agelength,wagepars,2)/2
    F[6]<-paramvals["F3"]*params$hours*PermWage(23*agelength,wagepars,2)/2
  }
  F[7] <- paramvals["Fold"]*params$hours*PermWage(23*agelength,wagepars,2)
  if(ptwork_halfcost==FALSE) F[8] <- paramvals["Fpartold"]*params$hours*PermWage(23*agelength,wagepars,2)
    else F[8] <- paramvals["Fold"]*params$hours*PermWage(23*agelength,wagepars,2)/2
  F<-F[!is.na(F)]
  
  if(single_samecost==FALSE){
    F_single<-NULL
    F_single[1]<-paramvals["F1_single"]*params$hours*PermWage(23*agelength,wagepars,2)
    F_single[2]<-paramvals["F2_single"]*params$hours*PermWage(23*agelength,wagepars,2)
    F_single[3]<-paramvals["F3_single"]*params$hours*PermWage(23*agelength,wagepars,2)
    if(ptwork_halfcost==FALSE){
      F_single[4]<-paramvals["Fpart1_single"]*params$hours*PermWage(23*agelength,wagepars,2)
      F_single[5]<-paramvals["Fpart2_single"]*params$hours*PermWage(23*agelength,wagepars,2)
      F_single[6]<-paramvals["Fpart3_single"]*params$hours*PermWage(23*agelength,wagepars,2)
    }
    else{
      F_single[4]<-paramvals["F1_single"]*params$hours*PermWage(23*agelength,wagepars,2)/2
      F_single[5]<-paramvals["F2_single"]*params$hours*PermWage(23*agelength,wagepars,2)/2
      F_single[6]<-paramvals["F3_single"]*params$hours*PermWage(23*agelength,wagepars,2)/2
    }
    F_single[7] <- paramvals["Fold_single"]*params$hours*PermWage(23*agelength,wagepars,2)
    if(ptwork_halfcost==FALSE) F_single[8] <- paramvals["Fpartold_single"]*params$hours*PermWage(23*agelength,wagepars,2)
    else F_single[8] <- paramvals["Fold_single"]*params$hours*PermWage(23*agelength,wagepars,2)/2
    
    F_single<-F_single[!is.na(F_single)]
  }
  else{
    F_single <- F
  }
  
  spF<-NULL
  spF[1]<-paramvals["spF1"]*params$hours*PermWage(23*agelength,spwagepars,2)
  spF[2]<-paramvals["spF2"]*params$hours*PermWage(23*agelength,spwagepars,2)
  if(ptwork_halfcost==FALSE){
    spF[3]<-paramvals["spFpart1"]*params$hours*PermWage(23*agelength,spwagepars,2)
    spF[4]<-paramvals["spFpart2"]*params$hours*PermWage(23*agelength,spwagepars,2)
  }
  else{
    spF[3]<-paramvals["spF1"]*params$hours*PermWage(23*agelength,spwagepars,2)/2
    spF[4]<-paramvals["spF2"]*params$hours*PermWage(23*agelength,spwagepars,2)/2
  }
  spF[5] <- paramvals["spFold"]*params$hours*PermWage(23*agelength,spwagepars,2)
  if(ptwork_halfcost==FALSE) spF[6] <- paramvals["spFpartold"]*params$hours*PermWage(23*agelength,wagepars,2)
  else spF[6] <- paramvals["spFold"]*params$hours*PermWage(23*agelength,spwagepars,2)/2
  spF<-spF[!is.na(spF)]
  
  param_input<-paramvals[which(names(paramvals)=="theta_mod"):max(which(names(paramvals)%in%c("delta","delta_single","spdelta")))]
  if(ptwork_halfdis==TRUE){
    param_input['etapart']<-param_input['eta']/2
    param_input['etapart_single']<-param_input['eta_single']/2
    param_input['spetapart']<-param_input['speta']/2
  }
  if(single_samedis==TRUE){
    param_input['eta_single']<-param_input['eta']
    param_input['etapt_single']<-param_input['etapart']
  }
  if(spouse_samejobrisk==TRUE) param_input["spdelta"]<-param_input["delta"]
  if(single_samejobrisk==TRUE) param_input["delta_single"]<-param_input["delta"]
  if(health_halfdis==TRUE) param_input['theta_sev']=param_input['theta_mod']*2
  param_input<-c(param_input,F,F_single,spF,paramvals[grepl("pi1.|pi2.|pi0.",names(paramvals))])
  names(param_input)<-names(paramvals)[which(names(paramvals)=="theta_mod"):length(paramvals)]
  if(single_sameallowanceprob==TRUE) param_input[grepl("pi1.|pi2.|pi0.",names(param_input)) & 
                                                 grepl("single",names(param_input))]<- param_input[
                                                   grepl("pi1.|pi2.|pi0.",names(param_input)) & 
                                                     grepl("married",names(param_input))]
  if(work_noallowanceprob==TRUE) param_input[grepl("pi1.|pi2.|pi0.",names(param_input)) &
                                               grepl("work",names(param_input))] <- 0
  
  if(single_sameprefs == TRUE) {
    param_input[c("eta_single","delta_single","F2_single","F3_single","Fold_single")] <- param_input[c("eta","delta","F2","F3","Fold")]
    param_input["spmu"]<-0
  }
  if(health_simpleprefs == TRUE) {
    param_input[c("theta_mod","eta_mod","F2","F2_single")] <- param_input[c("theta_sev","eta_sev","F3","F3_single")]/2
    param_input[c("Fold","Fold_single","spFold")]<-0
  }
  return(list(params=params,param_input=param_input))
}

paramsetup_old<-function(paramvals,wagepars,spwagepars,SSgenerosity=1,FSgenerosity=1,wageshifter=rep(1,50),spwageshifter=rep(1,50)){
  
  params<-list(
    agelength=agelength,
    wagevar= wagepars$wagevar,  #wage STDEV -- 0.027 #wagevar is variance of wage period shocks
    spwagevar = spwagepars$wagevar
    , tau = 1 #tau is the proportional tax used to fund social programs
    #THE FOLLOWING ARE PRE-SPECIFIED
    #    , reassess =  c(1-(1-0.222)^(4/agelength),1-(1-0.083)^(4/agelength),1-(1-0.036)^(4/agelength)) # reassess is the DI reassessment rate, vector,1 param per health level
    , reassess =  c(1,1/3,1/7) # reassess is the DI reassessment rate, vector,1 param per health level
    , gamma = 1.5
    , beta = 0.9756^(1/agelength)
    , R = 1.016^(1/agelength)
    ,wagereg_cons = c(0,0,0)
    ,wagereg_fe = c(wagepars$wagereg[[1]]$cons,wagepars$wagereg[[2]]$cons,wagepars$wagereg[[3]]$cons)
    ,wagereg_age1 = c(wagepars$wagereg[[1]]$age,wagepars$wagereg[[2]]$age,wagepars$wagereg[[3]]$age)
    ,wagereg_age2 = c(wagepars$wagereg[[1]]$agesq,wagepars$wagereg[[2]]$agesq,wagepars$wagereg[[3]]$agesq)
    ,wagereg_h1 = c(wagepars$wagereg[[1]]$health[1],wagepars$wagereg[[2]]$health[1],wagepars$wagereg[[2]]$health[1])
    ,wagereg_h2 = c(wagepars$wagereg[[1]]$health[2],wagepars$wagereg[[2]]$health[2],wagepars$wagereg[[2]]$health[2])
    ,wagereg_h3 = c(wagepars$wagereg[[1]]$health[3],wagepars$wagereg[[2]]$health[3],wagepars$wagereg[[2]]$health[3])
    ,spwagereg_cons = c(0,0,0)
    ,spwagereg_fe = c(spwagepars$wagereg[[1]]$cons,spwagepars$wagereg[[2]]$cons,spwagepars$wagereg[[3]]$cons)
    ,spwagereg_age1 = c(spwagepars$wagereg[[1]]$age,spwagepars$wagereg[[2]]$age,spwagepars$wagereg[[3]]$age)
    ,spwagereg_age2 = c(spwagepars$wagereg[[1]]$agesq,spwagepars$wagereg[[2]]$agesq,spwagepars$wagereg[[3]]$agesq)
    ,spwagereg_h1 = c(spwagepars$wagereg[[1]]$health[1],spwagepars$wagereg[[2]]$health[1],spwagepars$wagereg[[2]]$health[1])
    ,spwagereg_h2 = c(spwagepars$wagereg[[1]]$health[2],spwagepars$wagereg[[2]]$health[2],spwagepars$wagereg[[2]]$health[2]),
    hours=2000*agelength,
    kink1=5244/agelength, #formula cutoffs from 1996
    kink2=31620/agelength,
    kink3=62700/agelength,
    fsgross = 23920/agelength,
    fsnet= 18400/agelength, 
    UIcap = 12712/agelength,
    SSIval=5100/agelength,
    SSgenerosity=SSgenerosity,
    FSgenerosity=FSgenerosity,
    wageshifter=wageshifter,
    spwageshifter=spwageshifter
  )
  
  F<-NULL
  F[1]<-paramvals[6]*params$hours*PermWage(23*agelength,wagepars,2)
  F[2]<-paramvals[7]*params$hours*PermWage(23*agelength,wagepars,2)
  F[3]<-paramvals[8]*params$hours*PermWage(23*agelength,wagepars,2)
  
  spF<-NULL
  spF[1]<-paramvals[9]*params$hours*PermWage(23*agelength,spwagepars,2)
  spF[2]<-paramvals[10]*params$hours*PermWage(23*agelength,spwagepars,2)
  
  param_input<-paramvals
  param_input[6]<-F[1]
  param_input[7]<-F[2]
  param_input[8]<-F[3]
  param_input[9]<-spF[1]
  param_input[10]<-spF[2]
  return(list(params=params,param_input=param_input))
}

paramsetup_workint<-function(paramvals,wagepars,spwagepars,SSgenerosity=1,FSgenerosity=1,wageshifter=rep(1,50),spwageshifter=rep(1,50)){
  
  params<-list(
    agelength=agelength,
    #square root: sqrtm(sqrtm(params$health.prob))
    wagevar= wagepars$wagevar  #wage STDEV -- 0.027 #wagevar is variance of wage period shocks
    , tau = 1 #tau is the proportional tax used to fund social programs
    #THE FOLLOWING ARE PRE-SPECIFIED
    #    , reassess =  c(1-(1-0.222)^(4/agelength),1-(1-0.083)^(4/agelength),1-(1-0.036)^(4/agelength)) # reassess is the DI reassessment rate, vector,1 param per health level
    , reassess =  c(1,1/3,1/7) # reassess is the DI reassessment rate, vector,1 param per health level
    , gamma = 1.5
    , beta = 0.9756^(1/agelength)
    , R = 1.016^(1/agelength)
    ,wagereg_cons = c(0,0,0)
    ,wagereg_fe = c(wagepars$wagereg[[1]]$cons,wagepars$wagereg[[2]]$cons,wagepars$wagereg[[3]]$cons)
    ,wagereg_age1 = c(wagepars$wagereg[[1]]$age,wagepars$wagereg[[2]]$age,wagepars$wagereg[[3]]$age)
    ,wagereg_age2 = c(wagepars$wagereg[[1]]$agesq,wagepars$wagereg[[2]]$agesq,wagepars$wagereg[[3]]$agesq)
    ,wagereg_h1 = c(wagepars$wagereg[[1]]$health[1],wagepars$wagereg[[2]]$health[1],wagepars$wagereg[[2]]$health[1])
    ,wagereg_h2 = c(wagepars$wagereg[[1]]$health[2],wagepars$wagereg[[2]]$health[2],wagepars$wagereg[[2]]$health[2])
    ,wagereg_h3 = c(wagepars$wagereg[[1]]$health[3],wagepars$wagereg[[2]]$health[3],wagepars$wagereg[[2]]$health[3])
    ,spwagereg_cons = c(0,0,0)
    ,spwagereg_fe = c(spwagepars$wagereg[[1]]$cons,spwagepars$wagereg[[2]]$cons,spwagepars$wagereg[[3]]$cons)
    ,spwagereg_age1 = c(spwagepars$wagereg[[1]]$age,spwagepars$wagereg[[2]]$age,spwagepars$wagereg[[3]]$age)
    ,spwagereg_age2 = c(spwagepars$wagereg[[1]]$agesq,spwagepars$wagereg[[2]]$agesq,spwagepars$wagereg[[3]]$agesq)
    ,spwagereg_h1 = c(spwagepars$wagereg[[1]]$health[1],spwagepars$wagereg[[2]]$health[1],spwagepars$wagereg[[2]]$health[1])
    ,spwagereg_h2 = c(spwagepars$wagereg[[1]]$health[2],spwagepars$wagereg[[2]]$health[2],spwagepars$wagereg[[2]]$health[2]),
    hours=2000*agelength,
    kink1=4812/agelength, #formula cutoffs from 1993
    kink2=29040/agelength,
    kink3=66552/agelength,
    fsgross = 23920/agelength,
    fsnet= 18400/agelength, 
    UIcap = 12712/agelength,
    SSIval=5100/agelength,
    SSgenerosity=SSgenerosity,
    FSgenerosity=FSgenerosity,
    wageshifter=wageshifter,
    spwageshifter=spwageshifter
  )
  
  F<-NULL
  F[1]<-paramvals[7]*params$hours*PermWage(23*agelength,wagepars,2)
  F[2]<-paramvals[8]*params$hours*PermWage(23*agelength,wagepars,2)
  F[3]<-paramvals[9]*params$hours*PermWage(23*agelength,wagepars,2)
  
  spF<-NULL
  spF[1]<-paramvals[10]*params$hours*PermWage(23*agelength,spwagepars,2)
  spF[2]<-paramvals[11]*params$hours*PermWage(23*agelength,spwagepars,2)
  
  param_input<-paramvals
  param_input[7]<-F[1]
  param_input[8]<-F[2]
  param_input[9]<-F[3]
  param_input[10]<-spF[1]
  param_input[11]<-spF[2]
  return(list(params=params,param_input=param_input))
}




paramsetup_healthfe<-function(paramvals,wagepars,spwagepars,SSgenerosity=1,FSgenerosity=1,wageshifter=rep(1,50),spwageshifter=rep(1,50)){
  
  params<-list(
    agelength=agelength,
    #square root: sqrtm(sqrtm(params$health.prob))
    wagevar= wagepars$wagevar  #wage STDEV -- 0.027 #wagevar is variance of wage period shocks
    , tau = 1 #tau is the proportional tax used to fund social programs
    #THE FOLLOWING ARE PRE-SPECIFIED
    #    , reassess =  c(1-(1-0.222)^(4/agelength),1-(1-0.083)^(4/agelength),1-(1-0.036)^(4/agelength)) # reassess is the DI reassessment rate, vector,1 param per health level
    , reassess =  c(1,1/3,1/7) # reassess is the DI reassessment rate, vector,1 param per health level
    , gamma = 1.5
    , beta = 0.9756^(1/agelength)
    , R = 1.016^(1/agelength)
    ,wagereg_cons = c(0,0,0)
    ,wagereg_fe = c(wagepars$wagereg[[1]]$cons,wagepars$wagereg[[2]]$cons,wagepars$wagereg[[3]]$cons)
    ,wagereg_age1 = c(wagepars$wagereg[[1]]$age,wagepars$wagereg[[2]]$age,wagepars$wagereg[[3]]$age)
    ,wagereg_age2 = c(wagepars$wagereg[[1]]$agesq,wagepars$wagereg[[2]]$agesq,wagepars$wagereg[[3]]$agesq)
    ,wagereg_h1 = c(wagepars$wagereg[[1]]$health[1],wagepars$wagereg[[2]]$health[1],wagepars$wagereg[[2]]$health[1])
    ,wagereg_h2 = c(wagepars$wagereg[[1]]$health[2],wagepars$wagereg[[2]]$health[2],wagepars$wagereg[[2]]$health[2])
    ,wagereg_h3 = c(wagepars$wagereg[[1]]$health[3],wagepars$wagereg[[2]]$health[3],wagepars$wagereg[[2]]$health[3])
    ,spwagereg_cons = c(0,0,0)
    ,spwagereg_fe = c(spwagepars$wagereg[[1]]$cons,spwagepars$wagereg[[2]]$cons,spwagepars$wagereg[[3]]$cons)
    ,spwagereg_age1 = c(spwagepars$wagereg[[1]]$age,spwagepars$wagereg[[2]]$age,spwagepars$wagereg[[3]]$age)
    ,spwagereg_age2 = c(spwagepars$wagereg[[1]]$agesq,spwagepars$wagereg[[2]]$agesq,spwagepars$wagereg[[3]]$agesq)
    ,spwagereg_h1 = c(spwagepars$wagereg[[1]]$health[1],spwagepars$wagereg[[2]]$health[1],spwagepars$wagereg[[2]]$health[1])
    ,spwagereg_h2 = c(spwagepars$wagereg[[1]]$health[2],spwagepars$wagereg[[2]]$health[2],spwagepars$wagereg[[2]]$health[2]),
    hours=2000*agelength,
    kink1=4812/agelength, #formula cutoffs from 1993
    kink2=29040/agelength,
    kink3=66552/agelength,
    fsgross = 23920/agelength,
    fsnet= 18400/agelength, 
    UIcap = 12712/agelength,
    SSIval=5100/agelength,
    SSgenerosity=SSgenerosity,
    FSgenerosity=FSgenerosity,
    wageshifter=wageshifter,
    spwageshifter=spwageshifter
  )
  
  F<-NULL
  F[1]<-paramvals[6]*params$hours*PermWage(23*agelength,wagepars,2)
  F[2]<-paramvals[7]*params$hours*PermWage(23*agelength,wagepars,2)
  F[3]<-paramvals[8]*params$hours*PermWage(23*agelength,wagepars,2)
  
  Fint<-NULL
  Fint[1]<-0
  Fint[2]<-paramvals[9]*params$hours*PermWage(23*agelength,wagepars,2)
  
  spF<-NULL
  spF[1]<-paramvals[10]*params$hours*PermWage(23*agelength,spwagepars,2)
  spF[2]<-paramvals[11]*params$hours*PermWage(23*agelength,spwagepars,2)
  
  spFint<-NULL
  spFint[1]<-0
  spFint[2]<-paramvals[12]*params$hours*PermWage(23*agelength,spwagepars,2)
  spFint[2]<-paramvals[13]*params$hours*PermWage(23*agelength,spwagepars,2)
  
  param_input<-paramvals
  param_input[6]<-F[1]
  param_input[7]<-F[2]
  param_input[8]<-F[3]
  param_input[9]<-Fint[2]
  param_input[10]<-spF[1]
  param_input[11]<-spF[2]
  param_input[12]<-spFint[2]
  param_input[13]<-spFint[3]
  
  return(list(params=params,param_input=param_input))
}


const_moments<-function(output,shocks.wagenoise, shocks.spwagenoise,onlysingle=0,onlymarried=0,onlyFT=0,noflow=FALSE,nocomp=TRUE, nostockDI=TRUE,
                        nocontflow=TRUE,
                        wage=FALSE,
                        spwage=FALSE,
                        wagevar=FALSE,
                        cons=FALSE,
                        worksing=FALSE,
                        workmar=FALSE,
                        spwork=FALSE,
                        stockDI=FALSE,
                        compDI=FALSE,
                        flowDI=FALSE,
                        contflowDI=FALSE,
                        noeventDI=TRUE,
                        nohealthyDI=FALSE,
                        nospmoments=FALSE,
                        noworkDImoments=TRUE,
                        consmoments_simple=TRUE,
                        workmoments_simple=TRUE,
                        DImoments_simple=TRUE,
                        onlyevent=FALSE
          ){
  
  if(noeventDI==FALSE|onlyevent==TRUE){
  output[,age_firstdis:=min(age + 9999*(Health==1 | select!=1)),by=.(i,Type)]
  output[,age_firstmod:=min(age + 9999*(Health!=2 | select!=1)),by=.(i,Type)]
  output[,age_firstsev:=min(age + 9999*(Health!=3 | select!=1)),by=.(i,Type)]
  #Note: only people with age_firstdis>62 are people who are never work-limited
    #in observed years (select==1).
  output[,firstage:=min(age + 9999*(select!=1)),by=.(i,Type)]
  output[firstage==age_firstdis,age_firstdis:=9999]
  output[,firstdis_mod:=max((age==age_firstdis) & Health==2),by=.(i,Type)]
  output[,firstdis_sev:=max((age==age_firstdis) & Health==3),by=.(i,Type)]
  output[,mar_predis:=min((1-single) + 
                            9999*(age!=age_firstdis-2 | select != 1)),by=.(i,Type)]
  }
  
  output[,age2:=(age^2)/100]
  output[,married:=single==0]
  output[,anyWork:=Workpart==1 ]
  output[Work==1,anyWork:=NA ]
  output[,Workpart_single:=Workpart==1&single==1]
  output[,Workpart_married:=Workpart==1&single==0]
  output[,Work_single:=Work==1&single==1]
  output[,Work_married:=Work==1&single==0]
  output[,H2:=(Health==2)]
  output[,H3:=(Health==3)]
  output[,H1_single:=(Health==1)*single]
  output[,H2_single:=(Health==2)*single]
  output[,H3_single:=(Health==3)*single]
  output[,H1_married:=(Health==1)*(single==0)]
  output[,H2_married:=(Health==2)*(single==0)]
  output[,H3_married:=(Health==3)*(single==0)]
  output[,yearDI_single:=yearDI==1&single==1]
  output[,yearDI_H1_single:=yearDI==1&Health==1&single==1]
  output[,yearDI_H2_single:=yearDI==1&Health==2&single==1]
  output[,yearDI_H3_single:=yearDI==1&Health==3&single==1]
  output[,yearDI_married:=yearDI==1&single==0]
  output[,yearDI_H1_married:=yearDI==1&Health==1&single==0]
  output[,yearDI_H2_married:=yearDI==1&Health==2&single==0]
  output[,yearDI_H3_married:=yearDI==1&Health==3&single==0]
  output[single==1,spHealth:=1]
  output[,spH1:=spHealth==1]
  output[,spH2:=spHealth==2]
  output[,age_single:=age*single]
  output[,age_married:=age*(single==0)]
  output[,Work_Health:=Work*Health]
  output[,Workpart_Health:=Workpart*Health]
  output[,Work_H2:=Work==1&Health==2]
  output[,Work_H3:=Work==1&Health==3]
  output[,Workpart_H2:=Workpart==1&Health==2]
  output[,Workpart_H3:=Workpart==1&Health==3]
  output[,spWork_spH2:=spWork==1&spHealth==2]
  output[,spWorkpart_spH2:=spWorkpart==1&spHealth==2]
  output[,Type1:=Type==1]
  output[,Type2:=Type==2]
  output[,Type3:=Type==3]
  
  #Wage moments:
  output<-output[order(i,time),]
  output[,lw:=log(Wage) + shocks.wagenoise]
  output[Work!=1,lw:=NA]
  output[,sp_lw:=log(spWage) + shocks.spwagenoise]
  output[spWork!=1 | married!=1,sp_lw:=NA]
  #Comparing simulated vs actual wage reg, without selection correction
  output[,ageold_H2:= ageold*H2]
  output[,ageold_H3:= ageold*H3]
  output[,ageold_spH2:= ageold*spH2]
  output[select==1,age_limit:=min(age + 9999*(Health ==1)),by=.(i,Type)]
  output[select==1,age_work:=min(age + 9999*(Work ==0)),by=.(i,Type)]
  output[,insample:=select==1 & age_work < age_limit & age_work <9999]
  
  if(nrow(output[insample==1&Work==1,])==0){
  wagereg<-NA
  spwagereg<-NA
  mom.wage<-NA
  mom.spwage<-NA
  var.wage<-NA
  } else{
  wagereg<-tryCatch(feols(lw ~ as.numeric(Type1) + as.numeric(Type2) + as.numeric(Type3) + as.numeric(H3) + as.numeric(H2) + age + age2 + as.numeric(ageold) + as.numeric(ageold_H3) + as.numeric(ageold_H2) + as.numeric(married) - 1, data=output[insample==1&Work==1,]),error = function(e) "error")
  mom.wage <- tryCatch(wagereg$coefficients[names(wagereg$coefficients)!="as.numeric(married)"],error = function(e) "error")
  
  
  spwagereg <- tryCatch(feols(sp_lw ~ as.numeric(Type1) + as.numeric(Type2) + as.numeric(Type3) + as.numeric(spH2) + age + age2 + as.numeric(ageold) + as.numeric(ageold_spH2) - 1, data=output[insample==1&spWork==1 & married==1,]),error = function(e) "error")
  mom.spwage<-tryCatch(spwagereg$coefficients,error = function(e) "error")
  
  
  
  output[insample==1&Work==1,lwres:=tryCatch(wagereg$residuals,error = function(e) "error")]
  output[insample==1&spWork==1&married==1,sp_lwres:=tryCatch(spwagereg$residuals,error = function(e) "error")]
  output[,lwreslag:=c(NA,tryCatch(lwres[-.N],error = function(e) "error")),by=i]
  output[,sp_lwreslag:=c(NA,tryCatch(sp_lwres[-.N],error = function(e) "error")),by=i]
  
  var.wage <- tryCatch(unlist(output[insample==1 & Work == 1,.(var(lwres,na.rm=TRUE),
                                                               cov(lwres,lwreslag,use="pairwise"))]),error = function(e) "error")
  var.wage <- c(var.wage,
                tryCatch(unlist(output[insample==1 & spWork == 1,.(var(sp_lwres,na.rm=TRUE),
                                                                   cov(sp_lwres,sp_lwreslag,use="pairwise"))]),error = function(e) "error"))
  var.wage <- c(var.wage,
                tryCatch(unlist(output[insample==1 & Work == 1 & spWork == 1,.(cov(lwres,sp_lwres,use="pairwise"))]),error = function(e) Inf))
  }
  #mom.cons<-tryCatch(lm(log(yearcons) ~ as.factor(Health) + Workpart*married + Work*married + spWorkpart + spWork + as.factor(spHealth) +as.factor(Health)*(yearDI)+poly(age,2,raw=TRUE),data=output[age<=62&selectC==1,])$coefficients[c(4,6,5,13,14,2,3,10,15,16,7,8,9)],error = function(e) "error")
  if(consmoments_simple==FALSE){ 
    mom.cons<-tryCatch(lm(log(yearcons) ~ H2_single + H3_single +  Work_single + Workpart_single +
                          yearDI_H1_single + yearDI_H2_single + yearDI_H3_single +
                          married +
                          H2_married + H3_married  + Work_married + Workpart_married + Work_H2 + Work_H3 + 
                          Workpart_H2 + Workpart_H3 + 
                          spH2  + spWork + spWorkpart  +
                          spWork_spH2 + spWorkpart_spH2 +
                        yearDI_H1_married + yearDI_H2_married + yearDI_H3_married +
                        poly(age_single,2,raw=TRUE) + poly(age_married,2,raw=TRUE) ,data=output[age<=62&selectC==1 & insample==1,])$coefficients[2:25],error = function(e) "error")
  if(onlysingle==1) mom.cons <- mom.cons[c(1:7,13,14,15,16)]
  if(onlymarried==1) mom.cons<-mom.cons[9:24]
  } else {
    mom.cons<-tryCatch(lm(log(yearcons) ~ Health +  
                            Work + Workpart +
                            yearDI + 
                            married +
                            Work_Health + 
                            Workpart_Health + 
                            spH2  + spWork + spWorkpart  +
                            spWork_spH2 + spWorkpart_spH2 +
                            poly(age,2,raw=TRUE) ,data=output[age<=62&selectC==1 & insample==1,])$coefficients,error = function(e) "error")
    mom.cons<-mom.cons[!names(mom.cons)%in%c("marriedTRUE","poly(age, 2, raw = TRUE)1","poly(age, 2, raw = TRUE)2","(Intercept)")]
    if(onlysingle==1) mom.cons <- mom.cons[names(mom.cons)%in%c("Health","Work","Workpart","Work_Health","Workpart_Health","yearDI")]
    if(onlymarried==1) mom.cons<-mom.cons[names(mom.cons)%in%c("Health","Work","Workpart","Work_Health","Workpart_Health",
                                                               "spH2","spWork","spWorkpart","spWork_spH2","spWorkpart_spH2","yearDI")]
  }
  if(onlyFT==1) mom.cons<-mom.cons[!grepl("Workpart",names(mom.cons))]
  if(onlyFT==2) mom.cons<-mom.cons[!grepl("spWorkpart",names(mom.cons))]
  if(nospmoments==1) mom.cons<-mom.cons[!grepl("sp",names(mom.cons))]

  output<-output[order(Health,ageold),]
  mom.workmar<-NULL
  if(workmoments_simple==TRUE){
    mom.workmar<-as.vector(as.matrix(output[age<=62&insample==1,.(Work=mean(Work)),by=.(Health,ageold)][,c('Health','ageold'):=NULL]))
    if((onlyFT==0 | onlyFT == 2)) mom.workpart.mar<-as.vector(as.matrix(output[age<=62&insample==1,.(Workpart=mean(Workpart)),by=.(Health,ageold)][,c('Health','ageold'):=NULL]))
    if((onlyFT==0 | onlyFT == 2)) mom.workmar <- c(mom.work.mar,mom.workpart.mar)
    mom.worksing<-NULL
  } else{
  if(onlymarried==1 | onlysingle==0) mom.work.mar<-as.vector(as.matrix(output[age<=62&insample==1&single==0,.(Work=mean(Work)),by=.(Health,ageold)][,c('Health','ageold'):=NULL]))
  if((onlymarried==1 | onlysingle==0) & (onlyFT==0 | onlyFT == 2)) mom.workpart.mar<-as.vector(as.matrix(output[age<=62&insample==1&single==0,.(Workpart=mean(Workpart)),by=.(Health,ageold)][,c('Health','ageold'):=NULL]))
  if((onlymarried==1 | onlysingle==0) & (onlyFT==0 | onlyFT == 2)) mom.workmar <- c(mom.work.mar,mom.workpart.mar)
  if((onlymarried==1 | onlysingle==0) & onlyFT==1) mom.workmar<-mom.work.mar
  
  mom.worksing<-NULL
  if(onlysingle==1 | onlymarried==0) mom.work.sing<-as.vector(as.matrix(output[age<=62&insample==1&single==1,.(Work=mean(Work)),by=.(Health,ageold)][,c('Health','ageold'):=NULL]))
  if((onlysingle==1 | onlymarried==0) & (onlyFT==0 | onlyFT == 2)) mom.workpart.sing<-as.vector(as.matrix(output[age<=62&insample==1&single==1,.(Workpart=mean(Workpart)),by=.(Health,ageold)][,c('Health','ageold'):=NULL]))
  if((onlysingle==1 | onlymarried==0) & (onlyFT==0 | onlyFT == 2))  mom.worksing <- c(mom.work.sing,mom.workpart.sing)
  if((onlysingle==1 | onlymarried==0) & onlyFT == 1) mom.worksing <- mom.work.sing
  }
  mom.spwork<-NULL
  mom.spworkpart<-NULL
  output<-output[order(spHealth,ageold),]
  if(nospmoments==0){
  if(onlymarried==1 | onlysingle==0) mom.spwork<-as.vector(as.matrix(output[age<=62&insample==1&single==0,.(spWork=mean(spWork)),by=.(spHealth,ageold)][,c('spHealth','ageold'):=NULL]))
  if((onlymarried==1 | onlysingle==0) & onlyFT==0) mom.spworkpart<-as.vector(as.matrix(output[age<=62&insample==1&single==0,.(spWorkpart=mean(spWorkpart)),by=.(spHealth,ageold)][,c('spHealth','ageold'):=NULL]))
  if((onlymarried==1 | onlysingle==0) & onlyFT==0) mom.spwork <- c(mom.spwork,mom.spworkpart)
  }
  
  
  output<-output[order(-single,ageold,Health),]
  if(DImoments_simple==TRUE){
    mom.flowDI<-as.vector(as.matrix(output[age<=62&insample==1,.(flowDI=mean(enterDI,na.rm=TRUE)),by=.(Health,ageold)][,c('Health','ageold'):=NULL]))
    if(nohealthyDI) mom.flowDI<- mom.flowDI[-c(1,4)]
    mom.contflowDI<-as.vector(as.matrix(output[age<=62&insample==1&LastDI==1,.(flowDI=mean(yearDI,na.rm=TRUE)),by=.(Health,ageold)][,c('Health','ageold'):=NULL]))
    if(nohealthyDI) mom.constflowDI<- mom.constflowDI[-c(1,4)]
  } else {
  if(onlymarried == 0 & onlysingle == 0) {
    mom.flowDI<-as.vector(as.matrix(output[age<=62&insample==1,.(flowDI=mean(enterDI,na.rm=TRUE)),by=.(single,Health,ageold)][,c('single','Health','ageold'):=NULL]))
    if(nohealthyDI) mom.flowDI<- mom.flowDI[-c(1,4,7,10)]
    mom.contflowDI<-as.vector(as.matrix(output[age<=62&insample==1&LastDI==1,.(flowDI=mean(yearDI,na.rm=TRUE)),by=.(single,Health,ageold)][,c('single','Health','ageold'):=NULL]))
    if(nohealthyDI) mom.constflowDI<- mom.constflowDI[-c(1,4,7,10)]
  }
  if(onlymarried == 1 & onlysingle == 0) {
    mom.flowDI<-as.vector(as.matrix(output[age<=62&insample==1 & single == 0,.(flowDI=mean(enterDI,na.rm=TRUE)),by=.(single,Health,ageold)][,c('single','Health','ageold'):=NULL]))
    if(nohealthyDI) mom.flowDI<- mom.flowDI[-c(1,4)]
    mom.contflowDI<-as.vector(as.matrix(output[age<=62&insample==1&LastDI==1&single==0,.(flowDI=mean(yearDI,na.rm=TRUE)),by=.(single,Health,ageold)][,c('single','Health','ageold'):=NULL]))
    if(nohealthyDI) mom.constflowDI<- mom.constflowDI[-c(1,4)]
  }
  if(onlymarried == 0 & onlysingle == 1) {
    mom.flowDI<-as.vector(as.matrix(output[age<=62&insample==1 & single == 1,.(flowDI=mean(enterDI,na.rm=TRUE)),by=.(single,Health,ageold)][,c('single','Health','ageold'):=NULL]))
    if(nohealthyDI) mom.flowDI<- mom.flowDI[-c(1,4)]
    mom.contflowDI<-as.vector(as.matrix(output[age<=62&insample==1&LastDI==1&single==1,.(flowDI=mean(yearDI,na.rm=TRUE)),by=.(single,Health,ageold)][,c('single','Health','ageold'):=NULL]))
    if(nohealthyDI) mom.constflowDI<- mom.constflowDI[-c(1,4)]
  }
  }
  mom.stockDI<-NULL
  mom.compDI<-NULL
  
  output<-output[order(ageold,decreasing=c(FALSE)),]
  if(DImoments_simple==TRUE){
    mom.compDI<-c(mom.compDI,
                  output[age<=62&yearDI==1&insample==1,.(DIcomp3=mean(Health==3)),by=.(ageold)]$DIcomp3,
                  output[age<=62&yearDI==1&insample==1,.(DIcomp2=mean(Health==2)),by=.(ageold)]$DIcomp2,
                  output[age<=62&yearDI==1&insample==1,.(DIcomp1=mean(Health==1)),by=.(ageold)]$DIcomp1 )
    if(nohealthyDI) mom.compDI<- mom.compDI[-c(5,6)]
  } else{
  if(onlysingle==1 | onlymarried==0) mom.compDI<-c(mom.compDI,
                                                   output[age<=62&yearDI==1&insample==1&single==1,.(DIcomp3=mean(Health==3)),by=.(ageold)]$DIcomp3,
                                                   output[age<=62&yearDI==1&insample==1&single==1,.(DIcomp2=mean(Health==2)),by=.(ageold)]$DIcomp2,
                                                   output[age<=62&yearDI==1&insample==1&single==1,.(DIcomp1=mean(Health==1)),by=.(ageold)]$DIcomp1 )
  if(onlymarried==1 | onlysingle==0) mom.compDI<-c(mom.compDI,
                                                   output[age<=62&yearDI==1&insample==1&single==0,.(DIcomp3=mean(Health==3)),by=.(ageold)]$DIcomp3,
                                                   output[age<=62&yearDI==1&insample==1&single==0,.(DIcomp2=mean(Health==2)),by=.(ageold)]$DIcomp2,
                                                   output[age<=62&yearDI==1&insample==1&single==0,.(DIcomp1=mean(Health==1)),by=.(ageold)]$DIcomp1 )
  if(nohealthyDI & onlysingle==0& onlymarried==0) mom.compDI<- mom.compDI[-c(5,6,11,12)]
  if(nohealthyDI & (onlysingle!=0| onlymarried!=0)) mom.compDI<- mom.compDI[-c(5,6)]
  }
  
  if(noworkDImoments==TRUE){
  output<-output[order(ageold,Health),]
  if(DImoments_simple==TRUE) {
    mom.stockDI<-c(mom.stockDI,as.vector(as.matrix(output[age<=62&insample==1,.(stockDI=mean(yearDI)),by=.(Health,ageold)][,c('Health','ageold'):=NULL])))
  } else{
  if(onlysingle==1 | onlymarried==0) mom.stockDI<-c(mom.stockDI,as.vector(as.matrix(output[age<=62&insample==1&single==1,.(stockDI=mean(yearDI)),by=.(Health,ageold)][,c('Health','ageold'):=NULL])))
  if(onlymarried==1 | onlysingle==0) mom.stockDI<-c(mom.stockDI,as.vector(as.matrix(output[age<=62&insample==1&single==0,.(stockDI=mean(yearDI)),by=.(Health,ageold)][,c('Health','ageold'):=NULL])))
    }
  if(nohealthyDI & onlysingle==0& onlymarried==0 & DImoments_simple!=TRUE) mom.stockDI<- mom.stockDI[-c(1,4,7,10)]
  if(nohealthyDI & (onlysingle!=0| onlymarried!=0 | DImoments_simple==TRUE)) mom.stockDI<- mom.stockDI[-c(1,4)]
  }
  else{
    output<-output[order(anyWork,ageold,Health),]
    if(DImoments_simple==TRUE) {
      mom.stockDI<-c(mom.stockDI,as.vector(as.matrix(output[age<=62&insample==1&!is.na(anyWork),.(stockDI=mean(yearDI*(anyWork==0))),by=.(Health,ageold)][,c('Health','ageold'):=NULL])))
      mom.stockDI<-c(mom.stockDI,as.vector(as.matrix(output[age<=62&insample==1&!is.na(anyWork),.(stockDI=mean(yearDI*(anyWork==1))),by=.(Health,ageold)][,c('Health','ageold'):=NULL])))
    } else{
    if(onlysingle==1 | onlymarried==0) mom.stockDI<-c(mom.stockDI,as.vector(as.matrix(output[age<=62&insample==1&single==1&!is.na(anyWork),.(stockDI=mean(yearDI*(anyWork==0))),by=.(Health,ageold)][,c('Health','ageold'):=NULL])))
    if(onlymarried==1 | onlysingle==0) mom.stockDI<-c(mom.stockDI,as.vector(as.matrix(output[age<=62&insample==1&single==0&!is.na(anyWork),.(stockDI=mean(yearDI*(anyWork==0))),by=.(Health,ageold)][,c('Health','ageold'):=NULL])))
    if(onlysingle==1 | onlymarried==0) mom.stockDI<-c(mom.stockDI,as.vector(as.matrix(output[age<=62&insample==1&single==1&!is.na(anyWork),.(stockDI=mean(yearDI*(anyWork==1))),by=.(Health,ageold)][,c('Health','ageold'):=NULL])))
    if(onlymarried==1 | onlysingle==0) mom.stockDI<-c(mom.stockDI,as.vector(as.matrix(output[age<=62&insample==1&single==0&!is.na(anyWork),.(stockDI=mean(yearDI*(anyWork==1))),by=.(Health,ageold)][,c('Health','ageold'):=NULL])))
    }
    if(nohealthyDI & onlysingle==0& onlymarried==0 & DImoments_simple!=TRUE) mom.stockDI<- mom.stockDI[-c(1,4,7,10,13,16,19,22)]
    if(nohealthyDI & (onlysingle!=0| onlymarried!=0| DImoments_simple==TRUE)) mom.stockDI<- mom.stockDI[-c(1,4,7,10)]
    }
  
  if(length(mom.cons)!=0) names(mom.cons)<-"cons"
  if((onlymarried==1 | onlysingle==0)&length(mom.workmar)!=0) names(mom.workmar)<-"work_mar"
  if((onlysingle==1 | onlymarried==0)&length(mom.worksing)!=0) names(mom.worksing)<-"work_sing"
  if((onlymarried==1 | onlysingle==0)&length(mom.spwork)!=0) names(mom.spwork)<-"spwork"
  if(length(mom.stockDI)!=0) names(mom.stockDI)<-"stockDI"
  if(length(mom.compDI)!=0) names(mom.compDI)<-"compDI"
  if(length(mom.flowDI)!=0) names(mom.flowDI)<-"flowDI"
  if(noflow==TRUE) mom.flowDI<-NULL
  if(nocomp==TRUE) mom.compDI<-NULL
  if(nostockDI==TRUE) mom.stockDI<-NULL
  if(nocontflow==TRUE) mom.contflowDI <- NULL
  
  if(onlyevent==TRUE){
    output[,event_time:=age-age_firstdis]
    dynamicevent<-output[age<=62&insample==1 
                         & age > age_firstdis - 5 & age < age_firstdis+5
                         & age_firstsev < age_firstdis+5,mean(yearDI),by=.(event_time,single)][order(single,event_time),]
  }
  if(noeventDI==FALSE){
    mom.eventDI<-NULL
    if(DImoments_simple==TRUE){
      mom.eventDI<-c(mom.eventDI,
                     output[age<=62&insample==1 
                            & age >= age_firstdis & age < age_firstdis+5
                            & age_firstsev >= age_firstdis+5,mean(yearDI)],
                     output[age<=62&insample==1 
                            & age >= age_firstdis & age < age_firstdis+5
                            & age_firstsev < age_firstdis+5,mean(yearDI)])
    } else{
    if(onlymarried==FALSE){
  mom.eventDI<-c(mom.eventDI,
                 output[age<=62&insample==1& mar_predis==0 
                  & age >= age_firstdis & age < age_firstdis+5
                  & age_firstsev >= age_firstdis+5,mean(yearDI)],
                  output[age<=62&insample==1& mar_predis==0 
                  & age >= age_firstdis & age < age_firstdis+5
                  & age_firstsev < age_firstdis+5,mean(yearDI)])
    }
    if(onlysingle==FALSE){
    mom.eventDI<-c(mom.eventDI, 
                   output[age<=62&insample==1& mar_predis==1 
                 & age >= age_firstdis & age < age_firstdis+5
                 & age_firstsev >= age_firstdis+5,mean(yearDI)],
                  output[age<=62&insample==1& mar_predis==1 
                & age >= age_firstdis & age < age_firstdis+5
                & age_firstsev < age_firstdis+5,mean(yearDI)])
    }
    }
  }
  else mom.eventDI<-NULL
  
  if(wage) moments<-c(mom.wage)
  else if(spwage) moments<-c(mom.spwage)
  else if(wagevar) moments<-c(var.wage)
  else if(cons) moments<-mom.cons
  else if(worksing) moments <- mom.worksing
  else if(workmar) moments <-mom.workmar
  else if(spwork) moments <- mom.spwork
  else if(stockDI) moments <- mom.stockDI
  else if(compDI) moments <- mom.compDI
  else if(flowDI) moments <- mom.flowDI
  else if(contflowDI) moments <- mom.contflowDI
  else if((onlysingle==0 & onlymarried==0) | workmoments_simple==TRUE)  moments<-c(mom.wage,
                                                      mom.spwage,
                                                      var.wage,
                                                      mom.cons,
                                                      mom.workmar,
                                                      mom.worksing,
                                                      mom.spwork,
                                                      mom.stockDI,
                                                      mom.compDI,
                                                      mom.flowDI,
                                                      mom.contflowDI,
                                                      mom.eventDI)
  else if(onlysingle==1) moments<-c(mom.wage,
                                    var.wage,
                                    mom.cons,
                                    mom.worksing,
                                    mom.stockDI,
                                    mom.compDI,
                                    mom.flowDI,
                                    mom.contflowDI,
                                    mom.eventDI)
  else if(onlymarried==1) moments<-c(mom.wage,
                                     mom.spwage,
                                     var.wage,
                                     mom.cons,
                                     mom.workmar,
                                     mom.spwork,
                                     mom.stockDI,
                                     mom.compDI,
                                     mom.flowDI,
                                     mom.contflowDI,
                                     mom.eventDI)
  if(onlyevent==TRUE) moments<-dynamicevent
  return(moments)
}




welfareC<- function(paramvals=c(theta=c(-0.448,-0.448*2), sptheta=-0.448,eta=c(-0.185), speta=c(-0.185),
                                delta=0.062
                                ,F1=0,F2=0.547,F3=0.952,spF1=0.25, spF2=0.952,
                               pi0.young=0.006,pi1.young=0.171,pi2.young=0.331,pi0.old=0.075,
                               pi1.old=0.180,pi2.old=0.626),
                    parnames,
                    particle=1,
                    ptwork_halfcost=FALSE,
                    single_samecost=FALSE,
                    ptwork_halfdis=FALSE,
                    single_samedis=FALSE,
                    spouse_samejobrisk=FALSE,
                    single_sameallowanceprob=FALSE,
                    work_noallowanceprob=FALSE,
                    single_sameprefs = TRUE,
                    health_simpleprefs = TRUE,
                    consmoments_simple = TRUE,
                    workmoments_simple = TRUE,
                    DImoments_simple = TRUE,
                    outfile='outputvals',specification,datamoments,weightmat,healthprob,marriageprob,foodstamps,
                              numkids,Asset,Wage, spWage,agelength,
                    wagenoise, spwagenoise,
                    shocks.wage, shocks.spwage, shocks.health, shocks.allow, shocks.reassess,
                    shocks.jobloss, shocks.spjobloss, shocks.marriage, marriageprob_init, numsims, lumpsum=0, select, selectC,seed,onlysingle,onlymarried,
                    onlyFT, workDI=0, noSSDI=0, noflow=FALSE,nocomp=FALSE, nostockDI = FALSE,nocontflow=FALSE,
                    noeventDI=FALSE,
                    nospmoments=FALSE, noworkDImoments=FALSE){
  
  
  temp<-wagesetup_fit(paramvals,0,0,wage_gridsize=10,spwage_gridsize=5#,filesuffix=typename
                      )
  Wage<-temp$Wage
  spWage<-temp$spWage
  wagepars<-temp$wagepars
  spwagepars<-temp$spwagepars
  shocks.wagereal <- exp(shocks.wage*wagepars$wagevar)
  shocks.spwagereal <- exp((shocks.spwage*(1-wagepars$wagecorr^2)^0.5 + shocks.wage*wagepars$wagecorr)*spwagepars$wagevar)
  
  shocks.wagenoise <- wagenoise*wagepars$noisevar
  shocks.spwagenoise <- spwagenoise*spwagepars$noisevar
  #spwagepars$wagevar<-0
  
  rm(temp)
    names(paramvals) = parnames
  if(specification=='base'){
    temp<-paramsetup(paramvals,wagepars,spwagepars,
                     ptwork_halfcost=ptwork_halfcost,
                     single_samecost=single_samecost,
                     ptwork_halfdis=ptwork_halfdis,
                     single_samedis=single_samedis,
                     spouse_samejobrisk=spouse_samejobrisk,
                     single_sameprefs = single_sameprefs,
                     health_simpleprefs = health_simpleprefs,
                     single_sameallowanceprob=single_sameallowanceprob,
                     work_noallowanceprob=work_noallowanceprob)
  }  
  if(specification=='healthfe'){
  temp<-paramsetup_healthfe(paramvals,wagepars,spwagepars)
  }  
  if(specification=='workint'){
  temp<-paramsetup_workint(paramvals,wagepars,spwagepars)
  }
  params<-temp$params
  param_input<-temp$param_input
  rm(temp)
  simdata<-Valsim_solve(as.vector(param_input),params,healthprob, marriageprob,foodstamps,numkids,oecd,aperm(Asset),aperm(Wage), aperm(spWage),
                        shocks.wagereal,shocks.spwagereal,shocks.health,shocks.allow,shocks.reassess,shocks.jobloss,shocks.spjobloss, shocks.marriage,
                        marriageprob_init,numsims,lumpsum=lumpsum,print=0,
                        onlysingle=onlysingle,onlymarried=onlymarried,onlyFT=onlyFT, workDI=workDI,noSSDI=noSSDI)
  simdata<-as.data.table(simdata)
  simdata[,Health:=Health+1]
  simdata[,spHealth:=spHealth+1]
  
  print(simdata[C<=0])
  simdata[,row:=1:nrow(simdata)]
  simdata[,Ceq:=C/oecd[[Type+1]][time+1,single+1],by=row]
  simdata[,age:=floor(time/params$agelength+23)]
  simdata[,ageold:=age>=45*params$agelength]
  simdata[,yearcons:=sum(Ceq),by=.(age,i,Type)]
  simdata[,yearDI:=DI==2 & age<=62*params$agelength]
  
  simdata[,select:=select]
  simdata[,selectC:=selectC]
  simdata[select==1,enterDI:=yearDI==1&c(NA,yearDI[-.N])==0&time<62*params$agelength,by=i]
  simdata[select==1,LastDI:=c(NA,yearDI[-.N])==1&time<62*params$agelength,by=i]
  simdata[,ApplyLast:=c(NA,Apply[-.N]),by=i]
  simdata[,RejectedLast:=c(0,Rejected[-.N]),by=i]
  simdata[,HealthLast:=c(NA,Health[-.N]),by=i]
  simdata[select==1&yearDI==1&enterDI!=1,enterDI:=NA]
  simdata[,Type:=Type+1]
  #  simdata[,transfers:=transfers(time,DI==2,age>62,Work*Wage*params$hours+spWork*spWage,Work*Wage*params$hours/exp(wagepars$wagereg[[1]]$health[Health]),params,foodstamps,wagepars,spwagepars,Type),by=.(row)]
  
  
  
  
#   solution<-Welfare_solve(as.vector(param_input),datamoments,healthprob,foodstamps,numkids,oecd,aperm(Asset),aperm(Wage),0,0)
#   
#   #PULL OUTPUT FROM SOLUTION:
#   
  # V<-aperm(array(unlist(solution$V),c(2,3,3,2,2,dim(Asset)[3],dim(Wage)[3],50*agelength,numtypes)))
  # C<-aperm(array(unlist(solution$C),c(2,3,3,2,2,dim(Asset)[3],dim(Wage)[3],50*agelength,numtypes)))
  # Savings<-aperm(array(unlist(solution$Savings),c(2,3,3,2,2,dim(Asset)[3],dim(Wage)[3],50*agelength,numtypes)))
  # Work<-aperm(array(unlist(solution$Work),c(2,3,3,2,2,dim(Asset)[3],dim(Wage)[3],50*agelength,numtypes)))
  # Apply<-aperm(array(unlist(solution$Apply),c(2,3,3,2,2,dim(Asset)[3],dim(Wage)[3],50*agelength,numtypes)))
  # 
  # V_work_noworksp<-aperm(array(unlist(solution$V_work_noworksp),c(2,3,dim(Asset)[3],dim(Wage)[3],50*agelength,numtypes)))
  # V_nowork_elig_noworksp<-aperm(array(unlist(solution$V_nowork_elig_noworksp),c(2,3,dim(Asset)[3],dim(Wage)[3],50*agelength,numtypes)))
  # V_nowork_inelig_noworksp<-aperm(array(unlist(solution$V_nowork_inelig_noworksp),c(2,3,dim(Asset)[3],dim(Wage)[3],50*agelength,numtypes)))
  # V_apply_noworksp<-aperm(array(unlist(solution$V_apply_noworksp),c(2,3,dim(Asset)[3],dim(Wage)[3],50*agelength,numtypes)))
  # V_DI_noworksp<-aperm(array(unlist(solution$V_DI_noworksp),c(2,3,dim(Asset)[3],dim(Wage)[3],50*agelength,numtypes)))
  # 
  # V_work_worksp<-aperm(array(unlist(solution$V_work_worksp),c(2,3,dim(Asset)[3],dim(Wage)[3],50*agelength,numtypes)))
  # V_nowork_elig_worksp<-aperm(array(unlist(solution$V_nowork_elig_worksp),c(2,3,dim(Asset)[3],dim(Wage)[3],50*agelength,numtypes)))
  # V_nowork_inelig_worksp<-aperm(array(unlist(solution$V_nowork_inelig_worksp),c(2,3,dim(Asset)[3],dim(Wage)[3],50*agelength,numtypes)))
  # V_apply_worksp<-aperm(array(unlist(solution$V_apply_worksp),c(2,3,dim(Asset)[3],dim(Wage)[3],50*agelength,numtypes)))
  # V_DI_worksp<-aperm(array(unlist(solution$V_DI_worksp),c(2,3,dim(Asset)[3],dim(Wage)[3],50*agelength,numtypes)))
  # 
  # 
  # C_work_noworksp<-aperm(array(unlist(solution$C_work_noworksp),c(2,3,dim(Asset)[3],dim(Wage)[3],50*agelength,numtypes)))
  # C_nowork_elig_noworksp<-aperm(array(unlist(solution$C_nowork_elig_noworksp),c(2,3,dim(Asset)[3],dim(Wage)[3],50*agelength,numtypes)))
  # C_nowork_inelig_noworksp<-aperm(array(unlist(solution$C_nowork_inelig_noworksp),c(2,3,dim(Asset)[3],dim(Wage)[3],50*agelength,numtypes)))
  # C_apply_noworksp<-aperm(array(unlist(solution$C_apply_noworksp),c(2,3,dim(Asset)[3],dim(Wage)[3],50*agelength,numtypes)))
  # C_DI_noworksp<-aperm(array(unlist(solution$C_DI_noworksp),c(2,3,dim(Asset)[3],dim(Wage)[3],50*agelength,numtypes)))
  # 
  # C_work_worksp<-aperm(array(unlist(solution$C_work_worksp),c(2,3,dim(Asset)[3],dim(Wage)[3],50*agelength,numtypes)))
  # C_nowork_elig_worksp<-aperm(array(unlist(solution$C_nowork_elig_worksp),c(2,3,dim(Asset)[3],dim(Wage)[3],50*agelength,numtypes)))
  # C_nowork_inelig_worksp<-aperm(array(unlist(solution$C_nowork_inelig_worksp),c(2,3,dim(Asset)[3],dim(Wage)[3],50*agelength,numtypes)))
  # C_apply_worksp<-aperm(array(unlist(solution$C_apply_worksp),c(2,3,dim(Asset)[3],dim(Wage)[3],50*agelength,numtypes)))
  # C_DI_worksp<-aperm(array(unlist(solution$C_DI_worksp),c(2,3,dim(Asset)[3],dim(Wage)[3],50*agelength,numtypes)))
#   
#   #####################
#   
#   #checking formally if V is monotonic in A:
#   for(typ in 1:dim(V)[1]){
#       for(t in 1:dim(V)[2]){
#       for(w in 1:dim(V)[3]){
#         for(emp in 1:dim(V)[5]){
#           for(spemp in 1:dim(V)[6]){
#             for(DI in 1:dim(V)[7]){  
#             for(h in 1:dim(V)[8]){  
#               for(sph in 1:dim(V)[9]){  
#                 if(all(V[typ,t,w,,emp,spemp,DI,h,sph]==cummax(V[typ,t,w,,emp,spemp,DI,h,sph]),na.rm=TRUE)==FALSE) print(paste('MONOTONICITY PROBLEM AT DIM',typ,t,w,emp,spemp,DI,h,sph))
#               }
#             }
#           }
#         }
#       }
#       }
#     }
#   }
#   
#   #checking formally if V is monotonic in W (outside retirement):
#   for(typ in 1:dim(V)[1]){
#       for(t in 1:dim(V)[2]){
#       for(A in 1:dim(V)[4]){
#         for(emp in 1:dim(V)[5]){
#           for(spemp in 1:dim(V)[6]){
#             for(DI in 1:dim(V)[7]){  
#             for(h in 1:dim(V)[8]){  
#               for(sph in 1:dim(V)[9]){  
#                 if(all(round(V[typ,t,,A,emp,spemp,DI,h,sph],5)==cummax(round(V[typ,t,,A,emp,spemp,DI,h,sph],5)),na.rm=TRUE)==FALSE) print(paste('MONOTONICITY PROBLEM AT DIM',typ,t,A,emp,spemp,DI,h,sph))
#               }
#             }
#           }
#         }
#       }
#     }
#     }
#   }
#   
#   
#   
#   #Same for choice-specific value functions:  
#   for(typ in 1:dim(V_nowork_inelig_noworksp)[1]){
#       for(t in 1:dim(V_nowork_inelig_noworksp)[2]){
#       for(w in 1:dim(V_nowork_inelig_noworksp)[3]){
#         for(h in 1:dim(V_nowork_inelig_noworksp)[5]){  
#           for(sph in 1:dim(V_nowork_inelig_noworksp)[6]){  
#           if(all(V_work_noworksp[typ,t,w,,h,sph]==cummax(V_work_noworksp[typ,t,w,,h,sph]),na.rm=TRUE)==FALSE) print(paste('MONOTONICITY PROBLEM, WORK AT DIM',typ,t,w,h,sph))
#           if(all(V_nowork_elig_noworksp[typ,t,w,,h,sph]==cummax(V_nowork_elig_noworksp[typ,t,w,,h,sph]),na.rm=TRUE)==FALSE) print(paste('MONOTONICITY PROBLEM, NO WORK AT DIM',typ,t,w,h,sph))
#           if(all(V_nowork_inelig_noworksp[typ,t,w,,h,sph]==cummax(V_nowork_inelig_noworksp[typ,t,w,,h,sph]),na.rm=TRUE)==FALSE) print(paste('MONOTONICITY PROBLEM, NO WORK AT DIM',typ,t,w,h,sph))
#           if(all(V_apply_noworksp[typ,t,w,,h,sph]==cummax(V_apply_noworksp[typ,t,w,,h,sph]),na.rm=TRUE)==FALSE) print(paste('MONOTONICITY PROBLEM, APPLY AT DIM',typ,t,w,h,sph))
#           if(all(V_DI_noworksp[typ,t,w,,h,sph]==cummax(V_DI_noworksp[typ,t,w,,h,sph]),na.rm=TRUE)==FALSE) print(paste('MONOTONICITY PROBLEM, DI AT DIM',typ,t,w,h,sph))
#             
#             if(all(V_work_worksp[typ,t,w,,h,sph]==cummax(V_work_worksp[typ,t,w,,h,sph]),na.rm=TRUE)==FALSE) print(paste('MONOTONICITY PROBLEM, WORK AT DIM',typ,t,w,h,sph))
#             if(all(V_nowork_elig_worksp[typ,t,w,,h,sph]==cummax(V_nowork_elig_worksp[typ,t,w,,h,sph]),na.rm=TRUE)==FALSE) print(paste('MONOTONICITY PROBLEM, NO WORK AT DIM',typ,t,w,h,sph))
#             if(all(V_nowork_inelig_worksp[typ,t,w,,h,sph]==cummax(V_nowork_inelig_worksp[typ,t,w,,h,sph]),na.rm=TRUE)==FALSE) print(paste('MONOTONICITY PROBLEM, NO WORK AT DIM',typ,t,w,h,sph))
#             if(all(V_apply_worksp[typ,t,w,,h,sph]==cummax(V_apply_worksp[typ,t,w,,h,sph]),na.rm=TRUE)==FALSE) print(paste('MONOTONICITY PROBLEM, APPLY AT DIM',typ,t,w,h,sph))
#             if(all(V_DI_worksp[typ,t,w,,h,sph]==cummax(V_DI_worksp[typ,t,w,,h,sph]),na.rm=TRUE)==FALSE) print(paste('MONOTONICITY PROBLEM, DI AT DIM',typ,t,w,h,sph))
#             
#           }
#         }
#       }
#     }
#   }
#   
#   #checking formally if V is monotonic in W (outside retirement):
#   for(typ in 1:dim(V_nowork_inelig_noworksp)[1]){
#       for(t in 1:dim(V_nowork_inelig_noworksp)[2]){
#       for(A in 1:dim(V_nowork_inelig_noworksp)[4]){
#         for(h in 1:dim(V_nowork_inelig_noworksp)[5]){  
#           for(sph in 1:dim(V_nowork_inelig_noworksp)[6]){  
#             if(all(round(V_work_noworksp[typ,t,,A,h,sph],5)==cummax(round(V_work_noworksp[typ,t,,A,h,sph],5)),na.rm=TRUE)==FALSE) print(paste('MONOTONICITY PROBLEM, WORK AT DIM',typ,t,A,h,sph))
#           if(all(round(V_nowork_inelig_noworksp[typ,t,,A,h,sph],5)==cummax(round(V_nowork_inelig_noworksp[typ,t,,A,h,sph],5)),na.rm=TRUE)==FALSE) print(paste('MONOTONICITY PROBLEM, NO WORK AT DIM',typ,t,A,h,sph))
#           if(all(round(V_nowork_elig_noworksp[typ,t,,A,h,sph],5)==cummax(round(V_nowork_elig_noworksp[typ,t,,A,h,sph],5)),na.rm=TRUE)==FALSE) print(paste('MONOTONICITY PROBLEM, NO WORK AT DIM',typ,t,A,h,sph))
#           if(all(round(V_apply_noworksp[typ,t,,A,h,sph],5)==cummax(round(V_apply_noworksp[typ,t,,A,h,sph],5)),na.rm=TRUE)==FALSE) print(paste('MONOTONICITY PROBLEM, APPLY AT DIM',typ,t,A,h,sph))
#           if(all(round(V_DI_noworksp[typ,t,,A,h,sph],5)==cummax(round(V_DI_noworksp[typ,t,,A,h,sph],5)),na.rm=TRUE)==FALSE) print(paste('MONOTONICITY PROBLEM, DI AT DIM',typ,t,A,h,sph))
#             
#             if(all(round(V_work_worksp[typ,t,,A,h,sph],5)==cummax(round(V_work_worksp[typ,t,,A,h,sph],5)),na.rm=TRUE)==FALSE) print(paste('MONOTONICITY PROBLEM, WORK AT DIM',typ,t,A,h,sph))
#             if(all(round(V_nowork_inelig_worksp[typ,t,,A,h,sph],5)==cummax(round(V_nowork_inelig_worksp[typ,t,,A,h,sph],5)),na.rm=TRUE)==FALSE) print(paste('MONOTONICITY PROBLEM, NO WORK AT DIM',typ,t,A,h,sph))
#             if(all(round(V_nowork_elig_worksp[typ,t,,A,h,sph],5)==cummax(round(V_nowork_elig_worksp[typ,t,,A,h,sph],5)),na.rm=TRUE)==FALSE) print(paste('MONOTONICITY PROBLEM, NO WORK AT DIM',typ,t,A,h,sph))
#             if(all(round(V_apply_worksp[typ,t,,A,h,sph],5)==cummax(round(V_apply_worksp[typ,t,,A,h,sph],5)),na.rm=TRUE)==FALSE) print(paste('MONOTONICITY PROBLEM, APPLY AT DIM',typ,t,A,h,sph))
#             if(all(round(V_DI_worksp[typ,t,,A,h,sph],5)==cummax(round(V_DI_worksp[typ,t,,A,h,sph],5)),na.rm=TRUE)==FALSE) print(paste('MONOTONICITY PROBLEM, DI AT DIM',typ,t,A,h,sph))
#             
#         }
#       }
#     }
#       }
#   }
#   print("Simulating...")
# 
#   Ntypes<-c(800,800)
#   output<-data.table()
#   
#   # ex<- Filter(function(x) is.function(get(x,.GlobalEnv)),ls(.GlobalEnv))
#   # clusterExport(cl,ex)
#   # output<-foreach(sim = 1:(Ntypes[1]+Ntypes[2]+Ntypes[3]),.combine='rbind',
#   #                 .noexport = c("interpolate2d")){
# 
# 
#   for(sim in 1:(Ntypes[1]+Ntypes[2])){
#     if(sim <=Ntypes[1]) type <-1
#   if(sim >Ntypes[1]&sim<=Ntypes[2]+Ntypes[1]) type <-2
# #  seed<-as.integer(rnorm(1)*10000)
# 
#   
#   output<-rbind(output,data.table(i=sim,simulate.choices(V=V[type,,,,,,,,],C=C[type,,,,,,,,],
#                                                          V_work_worksp=V_work_worksp[type,,,,,],V_nowork_inelig_worksp=V_nowork_inelig_worksp[type,,,,,],
#                                                          V_nowork_elig_worksp=V_nowork_elig_worksp[type,,,,,],
#                                                          V_apply_worksp=V_apply_worksp[type,,,,,],V_DI_worksp=V_DI_worksp[type,,,,,],
#                                                          V_work_noworksp=V_work_noworksp[type,,,,,],V_nowork_inelig_noworksp=V_nowork_inelig_noworksp[type,,,,,],
#                                                          V_nowork_elig_noworksp=V_nowork_elig_noworksp[type,,,,,],
#                                                          V_apply_noworksp=V_apply_noworksp[type,,,,,],V_DI_noworksp=V_DI_noworksp[type,,,,,],
#                                                          C_work_worksp=C_work_worksp[type,,,,,],C_nowork_inelig_worksp=C_nowork_inelig_worksp[type,,,,,],
#                                                          C_nowork_elig_worksp=C_nowork_elig_worksp[type,,,,,],
#                                                          C_apply_worksp=C_apply_worksp[type,,,,,],C_DI_worksp=C_DI_worksp[type,,,,,],
#                                                          C_work_noworksp=C_work_noworksp[type,,,,,],C_nowork_inelig_noworksp=C_nowork_inelig_noworksp[type,,,,,],
#                                                          C_nowork_elig_noworksp=C_nowork_elig_noworksp[type,,,,,],
#                                                          C_apply_noworksp=C_apply_noworksp[type,,,,,],C_DI_noworksp=C_DI_noworksp[type,,,,,],
#                                                          Work=Work[type,,,,,,,],Apply=Apply[type,,,,,,,],Asset=Asset[type,,],Wage=Wage[type,,],
#                                                          shocks.wage=shocks.wage[1:50+(sim-1)*50], shocks.health=shocks.health[1:50+(sim-1)*50],
#                                                          shocks.allow=shocks.allow[1:50+(sim-1)*50], shocks.reassess=shocks.reassess[1:50+(sim-1)*50],
#                                                          shocks.jobloss=shocks.jobloss[1:50+(sim-1)*50], shocks.spjobloss=shocks.spjobloss[1:50+(sim-1)*50],
#                                                          params=params,transfers=transfers,seed=seed,type=type))
#   )
# }
#   print(output[C<=0,])
#   
#   output[,ageold:=time>45*params$agelength]
#   #Making OECD adjustment to consumption
#   output[,row:=1:nrow(output)]
#   output[,Ceq:=C/oecd[[Type]][time-23*params$agelength+1],by=row]
#   #Aggregating consumption to the year-level
#   output[,age:=floor(time/params$agelength)]
#   output[,yearcons:=sum(Ceq),by=.(age,i,Type)]
#   #Aggregating DI receipt at year-level
#   output[,yearDI:=DI==3]
#   #Keeping the observations falling on the start of the year, "survey" period:
# #  regdata<-output[(time-1)/params$agelength==age,.(C,Health,spHealth,yearDI,Work,spWork,time,ageold,age)]
#   
#   output[,enterDI:=yearDI==1&c(NA,yearDI[-.N])==0&time<62*params$agelength]
#   output[yearDI==1&enterDI!=1,enterDI:=NA]
#   
  ###MOMENTS:
  #ADD SPOUSAL HEALTH TO CONSUMPTION REG
  #ADD SPOUSAL WORK RATES BY SPOUSAL HEALTH
  
  
  #CONSUMPTION:
  #coefficients from a regression of log consumption on
  #health, health-DI, DI, work, and age squared
  #and spousal health
  
  #EMPLOYMENT RATES:
  #share employed by health and age
  #share spouses employed by health and age
  
  #DISABILITY INSURANCE:
  #share of DI  recipients coming from each  health and age group
  #Share on DI  by health type and age
  #Flow rates onto DI by health type and age group

#  cbind(const_moments(simdata), const_moments(output), const_moments(simdata)-const_moments(output))
  moments <- tryCatch(const_moments(simdata,
                           shocks.wagenoise=shocks.wagenoise, shocks.spwagenoise=shocks.spwagenoise,
                           onlysingle=onlysingle,onlymarried=onlymarried,onlyFT=onlyFT,noflow=noflow,nocomp=nocomp,noeventDI=noeventDI, nostockDI=nostockDI, nocontflow=nocontflow, nospmoments=nospmoments,noworkDImoments=noworkDImoments,
                           consmoments_simple=consmoments_simple,DImoments_simple=DImoments_simple,workmoments_simple=workmoments_simple),
  error=function(e) "Error")
  #DEFINE datamoments -- 30 elements
  if(length(datamoments)!=length(moments)|sum(is.na(moments))!=0) objective<-Inf
  else tryCatch(objective<-t(moments-datamoments)%*%weightmat%*%as.matrix(moments-datamoments),
                error=function(e) "Error")
  write_csv(x=as.data.frame(t(c(particle,paramvals,objective,999))),path=paste0(outfile,".csv"),append=TRUE,col_names=FALSE)
  print(c(objective," from particle: ",particle, " and params: ",paramvals))
  return(objective)
}


  Wageval_solve<-function(V0,type,t,person,wage,spwage,asset,DI,emp,spemp,health,sphealth,V1,Wage,spWage,Asset, policyname,
                          wtpDI=0, tol=1e-5){
#    print(paste0("solving policy ",policyname," for time ", t," of person ", person))
  if(wtpDI==1) DI1<-2
  else DI1<-DI
  
  interpV<-function(V,DI){
    if(dim(spWage)[3]>1){
    point<-(wageinterp*spwageinterp*assetinterp*V[type,t+1,wagepoint+1,spwagepoint+1,assetpoint+1,emp+1,spemp+1,DI+1,health+1,sphealth+1]
    + (1-wageinterp)*spwageinterp*assetinterp*V[type,t+1,wagepoint,spwagepoint+1,assetpoint+1,emp+1,spemp+1,DI+1,health+1,sphealth+1]
    + wageinterp*(1-spwageinterp)*assetinterp*V[type,t+1,wagepoint+1,spwagepoint,assetpoint+1,emp+1,spemp+1,DI+1,health+1,sphealth+1]
    + wageinterp*spwageinterp*(1-assetinterp)*V[type,t+1,wagepoint+1,spwagepoint+1,assetpoint,emp+1,spemp+1,DI+1,health+1,sphealth+1]
    + (1-wageinterp)*(1-spwageinterp)*assetinterp*V[type,t+1,wagepoint,spwagepoint,assetpoint+1,emp+1,spemp+1,DI+1,health+1,sphealth+1]
    + (1-wageinterp)*spwageinterp*(1-assetinterp)*V[type,t+1,wagepoint,spwagepoint+1,assetpoint,emp+1,spemp+1,DI+1,health+1,sphealth+1]
    + wageinterp*(1-spwageinterp)*(1-assetinterp)*V[type,t+1,wagepoint+1,spwagepoint,assetpoint,emp+1,spemp+1,DI+1,health+1,sphealth+1]
    + (1-wageinterp)*(1-spwageinterp)*(1-assetinterp)*V[type,t+1,wagepoint,spwagepoint,assetpoint,emp+1,spemp+1,DI+1,health+1,sphealth+1])
    }
    else{
      point<-(wageinterp*(1-spwageinterp)*assetinterp*V[type,t+1,wagepoint+1,spwagepoint,assetpoint+1,emp+1,spemp+1,DI+1,health+1,sphealth+1]
              + (1-wageinterp)*(1-spwageinterp)*assetinterp*V[type,t+1,wagepoint,spwagepoint,assetpoint+1,emp+1,spemp+1,DI+1,health+1,sphealth+1]
              + wageinterp*(1-spwageinterp)*(1-assetinterp)*V[type,t+1,wagepoint+1,spwagepoint,assetpoint,emp+1,spemp+1,DI+1,health+1,sphealth+1]
              + (1-wageinterp)*(1-spwageinterp)*(1-assetinterp)*V[type,t+1,wagepoint,spwagepoint,assetpoint,emp+1,spemp+1,DI+1,health+1,sphealth+1])
    }
    return(point)
  }
  newwage<-wage
  spnewwage<-spwage
  wagepoint=max(c(1,which(newwage>Wage[type,t+1,-dim(Wage)[3]])))
  if(dim(spWage)[3]>1) spwagepoint=max(c(1,which(spnewwage>spWage[type,t+1,-dim(spWage)[3]])))
  else spwagepoint = 1
  assetpoint=max(c(1,which(asset>Asset[type,t+1,-dim(Asset)[3]])))
  wageinterp= min((newwage-Wage[type,t+1,wagepoint])/(Wage[type,t+1,wagepoint+1]-Wage[type,t+1,wagepoint]),1)
  if(dim(spWage)[3]>1) spwageinterp= min((spnewwage-spWage[type,t+1,spwagepoint])/(spWage[type,t+1,spwagepoint+1]-spWage[type,t+1,spwagepoint]),1)
  else spwageinterp = 0
  assetinterp= min((asset-Asset[type,t+1,assetpoint])/(Asset[type,t+1,assetpoint+1]-Asset[type,t+1,assetpoint]),1)

  V<-V0
  V_base<- interpV(V,DI)
  
  V<-V1
  V_new<- interpV(V,DI1)
  lowermult <- 0
  uppermult <- 1
  mult<-1
  while((V_new < V_base)
        & (wagepoint+wageinterp<dim(Wage)[3] | 
           spwagepoint + spwageinterp < dim(spWage)[3])){
  uppermult <- uppermult*2
  mult<-mult*2
  newwage <- wage*mult
  spnewwage<-spwage*mult
  wagepoint<-max(c(1,which(newwage>Wage[type,t+1,-dim(Wage)[3]])))
  spwagepoint=max(c(1,which(spnewwage>spWage[type,t+1,-dim(spWage)[3]])))
  wageinterp<- min((newwage-Wage[type,t+1,wagepoint])/(Wage[type,t+1,wagepoint+1]-Wage[type,t+1,wagepoint]),1)
  if(dim(spWage)[3]>1) spwageinterp= min((spnewwage-spWage[type,t+1,spwagepoint])/(spWage[type,t+1,spwagepoint+1]-spWage[type,t+1,spwagepoint]),1)
  else spwageinterp=0
  V_new<- interpV(V,DI1)
  }
  while((abs(V_new-V_base) > tol*(abs(V_new)+abs(V_base)))
        & (newwage>min(Wage[type,t+1,]))
        & (spnewwage > min(spWage[type,t+1,]))
        & (wagepoint+wageinterp<dim(Wage)[3])
        & (spwagepoint+spwageinterp<dim(spWage)[3])){

   if(V_new > V_base){
     lowermult <- lowermult
     uppermult <- mult
   } 
    else{
      lowermult <- mult
      upperwage <- uppermult
    }
    mult <- (uppermult-lowermult)/2+lowermult
    newwage<-wage*mult
    spnewwage<-spwage*mult
    wagepoint<-max(c(1,which(newwage>Wage[type,t+1,-dim(Wage)[3]])))
    spwagepoint=max(c(1,which(spnewwage>spWage[type,t+1,-dim(spWage)[3]])))
    wageinterp<- min((newwage-Wage[type,t+1,wagepoint])/(Wage[type,t+1,wagepoint+1]-Wage[type,t+1,wagepoint]),1)
    if(dim(spWage)[3]>1) spwageinterp= min((spnewwage-spWage[type,t+1,spwagepoint])/(spWage[type,t+1,spwagepoint+1]-spWage[type,t+1,spwagepoint]),1)
    else spwageinterp=0
    V_new<- interpV(V,DI1)
  }
  if(abs(V_new-V_base) <= tol*(abs(V_new)+abs(V_base))) output <- -(newwage - wage)*100/wage
  else if(V_new > V_base) output<- Inf
  else if(V_new < V_base) output<- -Inf
  
  return (output)
  
}



Assetval_solve<-function(V0,type,t,person,wage,spwage,asset,DI,emp,spemp,health,sphealth,V1,Wage,spWage,Asset, policyname,
                        wtpDI=0,tol=1e-5,asset_tol=100){
#  print(paste0("solving policy ",policyname," for time ", t," of person ", person))
  if(wtpDI==1) DI1<-2
  else DI1<-DI

  interpV<-function(V,DI){
    if(dim(spWage)[3]>1){
      point<-(wageinterp*spwageinterp*assetinterp*V[type,t+1,wagepoint+1,spwagepoint+1,assetpoint+1,emp+1,spemp+1,DI+1,health,sphealth]
              + (1-wageinterp)*spwageinterp*assetinterp*V[type,t+1,wagepoint,spwagepoint+1,assetpoint+1,emp+1,spemp+1,DI+1,health,sphealth]
              + wageinterp*(1-spwageinterp)*assetinterp*V[type,t+1,wagepoint+1,spwagepoint,assetpoint+1,emp+1,spemp+1,DI+1,health,sphealth]
              + wageinterp*spwageinterp*(1-assetinterp)*V[type,t+1,wagepoint+1,spwagepoint+1,assetpoint,emp+1,spemp+1,DI+1,health,sphealth]
              + (1-wageinterp)*(1-spwageinterp)*assetinterp*V[type,t+1,wagepoint,spwagepoint,assetpoint+1,emp+1,spemp+1,DI+1,health,sphealth]
              + (1-wageinterp)*spwageinterp*(1-assetinterp)*V[type,t+1,wagepoint,spwagepoint+1,assetpoint,emp+1,spemp+1,DI+1,health,sphealth]
              + wageinterp*(1-spwageinterp)*(1-assetinterp)*V[type,t+1,wagepoint+1,spwagepoint,assetpoint,emp+1,spemp+1,DI+1,health,sphealth]
              + (1-wageinterp)*(1-spwageinterp)*(1-assetinterp)*V[type,t+1,wagepoint,spwagepoint,assetpoint,emp+1,spemp+1,DI+1,health,sphealth])
    }
    else{
      point<-(wageinterp*(1-spwageinterp)*assetinterp*V[type,t+1,wagepoint+1,spwagepoint,assetpoint+1,emp+1,spemp+1,DI+1,health,sphealth]
              + (1-wageinterp)*(1-spwageinterp)*assetinterp*V[type,t+1,wagepoint,spwagepoint,assetpoint+1,emp+1,spemp+1,DI+1,health,sphealth]
              + wageinterp*(1-spwageinterp)*(1-assetinterp)*V[type,t+1,wagepoint+1,spwagepoint,assetpoint,emp+1,spemp+1,DI+1,health,sphealth]
              + (1-wageinterp)*(1-spwageinterp)*(1-assetinterp)*V[type,t+1,wagepoint,spwagepoint,assetpoint,emp+1,spemp+1,DI+1,health,sphealth])
    }
    return(point)
  }
  
  
  newasset<-asset
  assetpoint=max(c(1,which(newasset>Asset[type,t+1,-dim(Asset)[3]])))
  wagepoint=max(c(1,which(wage>Wage[type,t+1,-dim(Wage)[3]])))
  if(dim(spWage)[3]>1) spwagepoint=max(c(1,which(spwage>spWage[type,t+1,-dim(spWage)[3]])))
  else spwagepoint = 1
  
  
  wageinterp= min((wage-Wage[type,t+1,wagepoint])/(Wage[type,t+1,wagepoint+1]-Wage[type,t+1,wagepoint]),1)
  if(dim(spWage)[3]>1) spwageinterp= min((spwage-spWage[type,t+1,spwagepoint])/(spWage[type,t+1,spwagepoint+1]-spWage[type,t+1,spwagepoint]),1)
  else spwageinterp = 0
  assetinterp= min((newasset-Asset[type,t+1,assetpoint])/(Asset[type,t+1,assetpoint+1]-Asset[type,t+1,assetpoint]),1)
  V<-V0
  V_base<- interpV(V,DI)
  
  V<-V1
  V_new<- interpV(V,DI1)
  
  lowermult <- 0
  uppermult <- 1
  mult<-1
  while(V_new < V_base 
        & (assetpoint+assetinterp<dim(Asset)[3])){
    uppermult <- uppermult*2
    mult<-mult*2
    newasset<-(asset+1000)*mult
    assetpoint=max(c(1,which(newasset>Asset[type,t+1,-dim(Asset)[3]])))
    
    assetinterp= min((newasset-Asset[type,t+1,assetpoint])/(Asset[type,t+1,assetpoint+1]-Asset[type,t+1,assetpoint]),1)
    V_new<- interpV(V,DI1)
  }
  while((abs(V_new-V_base) > tol*(abs(V_new)+abs(V_base))) 
        & (newasset>asset_tol)
        & (assetpoint+assetinterp<dim(Asset)[3])){
    
    if(V_new > V_base){
      lowermult <- lowermult
      uppermult <- mult
    } 
    else{ 
      lowermult <- mult
      upperwage <- uppermult
    }
    mult <- (uppermult-lowermult)/2+lowermult
    newasset<-(asset+1000)*mult
    assetpoint=max(c(1,which(newasset>Asset[type,t+1,-dim(Asset)[3]])))
    assetinterp= min((newasset-Asset[type,t+1,assetpoint])/(Asset[type,t+1,assetpoint+1]-Asset[type,t+1,assetpoint]),1)
    V_new<- interpV(V,DI1)
  }
  if(abs(V_new-V_base) <= tol*(abs(V_new)+abs(V_base))) output <- (asset - newasset)
  else if(V_new > V_base) output<- Inf
  else if(V_new < V_base) output<- -Inf
  return (as.double(output))
  
}

networkinc<-function(potentialwage,potentialwagesp,work,worksp,single,t,numkids,type){
  pretaxinc<-potentialwage*work+potentialwagesp*worksp
  taxableinc<-max(pretaxinc-6700*(1-single)-4000*single-2550*(1 + (1 - single) + (1 - single) * numkids[[type]][t+1,single+1]),0)
  if(single==0){
    if (taxableinc <= 40100) taxsim = max(taxableinc * 0.15, 0.0)
    if (taxableinc > 40100 & taxableinc <= 96900) taxsim = 40100 * 0.15 + (taxableinc - 40100) * 0.28
    if (taxableinc > 96900 & taxableinc <= 147700) taxsim = 40100 * 0.15 + (96900 - 40100) * 0.28 + (taxableinc - 96900) * 0.31
    if (taxableinc > 147700 & taxableinc <= 263750) taxsim = 40100 * 0.15 + (96900 - 40100) * 0.28 + (147700 - 96900) * 0.31 + (taxableinc - 147700) * 0.36
    if (taxableinc > 263750) taxsim = 40100 * 0.15 + (96900 - 40100) * 0.28 + (147700 - 96900) * 0.31 + (263750 - 147700) * 0.36 + (taxableinc - 263750) * 0.396
  }
  if(single==1){
    if (taxableinc <= 24000) taxsim = max(taxableinc * 0.15, 0.0)
    if (taxableinc > 24000 & taxableinc <= 58150) taxsim = 24000 * 0.15 + (taxableinc - 24000) * 0.28
    if (taxableinc > 58150 & taxableinc <= 121300) taxsim = 24000 * 0.15 + (58150 - 24000) * 0.28 + (taxableinc - 58150) * 0.31
    if (taxableinc > 121300 & taxableinc <= 263750) taxsim = 24000 * 0.15 + (58150 - 24000) * 0.28 + (121300 - 58150) * 0.31 + (taxableinc - 121300) * 0.36
    if (taxableinc > 263750) taxsim = 24000 * 0.15 + (58150 - 24000) * 0.28 + (121300 - 58150) * 0.31 + (263750 - 121300) * 0.36 + (taxableinc - 263750) * 0.396
  }
  #Child tax credit:
  taxsim <- max(taxsim-400*numkids[[type]][t+1,single+1])
  #FICA taxes:
  if (work * potentialwage < 62700) taxsim = taxsim + 0.0765 * work * potentialwage
  if (work * potentialwage >= 62700) taxsim = taxsim + 62700 * 0.062 + 0.0145 * work * potentialwage
  if (worksp * potentialwagesp < 62700) taxsim = taxsim + 0.0765 * worksp * potentialwagesp
  if (worksp * potentialwagesp >= 62700) taxsim = taxsim + 62700 * 0.062 + 0.0145 * worksp * potentialwagesp
  eitc0<-0
  eitc1<-0
  eitc<-0
  #EITC
  if (numkids[[type]][t+1, single+1] <= 1 | single == 1) {
    if (pretaxinc <= 4200) eitc0 = pretaxinc * 0.0765
    if (pretaxinc > 4200 && pretaxinc <= 5300) eitc0 = 4200 * 0.0765
    if (pretaxinc > 5300 && pretaxinc <= 9500) eitc0 = 4200 * 0.0765 - (pretaxinc - 5300) * 0.0765
    if (pretaxinc <= 6350) eitc1 = pretaxinc * 0.34
    if (pretaxinc > 6350 && pretaxinc <= 11650) eitc1 = 6350 * 0.34
    if (pretaxinc > 11650 && pretaxinc <= 25100) eitc1 = 6350 * 0.34 - (pretaxinc - 11650) * 0.1598
    if (pretaxinc <= 25100) eitc = eitc0 * (1 - numkids[[type]][t+1, single+1]) + eitc1 * (numkids[[type]][t+1, single+1])
    if (single == 1) eitc = eitc0
  }
  if (numkids[[type]][t+1, single+1] > 1 & single == 0) {
    if (pretaxinc <= 6350) eitc0 = pretaxinc * 0.34
    if (pretaxinc > 6350 & pretaxinc <= 11650) eitc0 = 6350 * 0.34
    if (pretaxinc > 11650 & pretaxinc <= 25100) eitc0 = 6350 * 0.34 - (pretaxinc - 11650) * 0.1598
    if (pretaxinc <= 8890) eitc1 = pretaxinc * 0.4
    if (pretaxinc > 8890 & pretaxinc <= 11650) eitc1 = 8890 * 0.4
    if (pretaxinc > 11650 & pretaxinc <= 28495) eitc1 = 8890 * 0.4 - (pretaxinc - 11650) * 0.2106
    if (pretaxinc <= 28495) eitc = eitc0 * max(2 - numkids[[type]][t+1, single+1], 0.0) +
                                   eitc1 * min(numkids[[type]][t+1, single+1] - 1, 1.0)
  }
  taxsim = taxsim - eitc;
  netinc<-pretaxinc - taxsim 
  return(netinc)
  }

moralhazard<-function(Ceq,Work,spWork,single,Health,spHealth,wage,spwage,time,ageold,DI,Apply,numkids,Type,paramvals){
if((DI==0 & Apply == 0)|time+23 > 62) wouldwork=NA
else{
  netinc_DI=networkinc(wage,spwage,Work,spWork,single,time,numkids,Type)
  netinc_work=networkinc(wage,spwage,1,spWork,single,time,numkids,Type)
  newworkinc<-netinc_work-netinc_DI
  flow_DI<-flow(Ceq,
                Work,spWork,single,Health,spHealth,paramvals)
  flow_Work<-flow(Ceq+newworkinc-paramvals[["F1"]]*(Health==1&single==0)
                  -paramvals[["F2"]]*(Health==2&single==0)
                  -paramvals[["F3"]]*(Health==3&single==0)
                  -paramvals[["Fold"]]*(ageold==1&single==0)
                  -paramvals[["F1_single"]]*(Health==1&single==1)
                  -paramvals[["F2_single"]]*(Health==2&single==1)
                  -paramvals[["F3_single"]]*(Health==3&single==1)
                  -paramvals[["Fold_single"]]*(ageold==1&single==1),
                  1,spWork,single,Health,spHealth,paramvals)
  wouldwork <- flow_DI < flow_Work
}
  return(wouldwork)
}
